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ABSTRACT 


The polycondensation stage of polyethylene terephthalate 
(PET) formation is assumed to include side reactions leading 
to the formation of diethylene glycol, vinyl end groups j acid 
end groups^ in addition to the usual polymerization of his 
(2- hydroxyethyl) terephthalate (BHET) in semi-hatch reactors* 

A flexible objective function has been proposed "with temperature 
as well as pressure as control variables, 0oii5)utations show 
that the pressure should be reduced to the lowest limit under 
all possible conditions. Consequently optimal temperature 
profiles in batch reactors.. are obtained for various lower limits 
of reactor pressures using the combination of first and second 
variation techniques. For the first variation technique, the 
vector iteration method of computation was used and the near 
optimal profile s© obtained was used as the initial guess for the 
second variation technique. 

The result of optimization shows that the lower limit of 
pressure and weighting parameters appearing in the objective 
function have profound effect on the optimum profiles,. For 
higher pressures, one obtains that a high temperature must be 
used initially which must be lowered later on to minimize the 
foimiation of side products. However for lower pressures, the 
temperature must be increased from a low value inltiallyj 
however for large times of polymerization, this must be lowered 
to minimize the formation of side products. It is thus seen that 
the optimum temperature profile exhibits a broad maximum, 



CHAPTER 1 
I HDRODUCTIOH 

Among the saturated polyesters, polyethylene terephthalate 

^ ?T 

(PBa),4CH2CH2 - 0 - C-<o>- 0 - 0 4^ ) is the most important 

polymer from commercial point of view^ It is also loiown hy 

various other names such as ’’Dacron” and "Terylene”. Eibre— 

1 ? 

grade PET is mainly used for making textiles * and its; 
contribution to the world synthetic fibre demand is around 

f 

forty percent* It has many other industrial applications, for 
example it can be used as a raw material for the production 
of a new transparent film for electrical insulation and also as 
a base for a new photdgraphic film. It can also be used as a 
moulding material. 

Method of PET Production; 

In Industries, the production of fibre-grade PET is 
carried out in four stages.^ These are: (1) Transesterifieation 
or direct esterification, (2) Prepolymerization, (3) Melt 
polycondensation, and (4) Solid state polycondensation. The 
different stages are as shown in ^’ig, 1. The monomer for the 
production of PET is bis (hydroxyethyl) Terephthalate (H3ET, 

HOCH 2 CH 2 OG ~<o^G 0GH2CH20H) . In the first stage, this 

monomer is synthesized either through a transesterifieation 

e»r a direct esterification reaction, using Dimethyl Terephthalate 
0 0 

(IM, and Ethylene Glycol (EG, OHGH 2 GH 2 OH) 
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as the raw materials. In transesterification, U® and BJ 
are reacted and the condensation product, methanol is conti- 
nuously removed, whereas in direct esterification, terephthalic 
acid (TPA) is employed in place of DI®, and the reaction takes 
place between this and with the evolution of water 

as the condensation product, which is removed from the reaction 
mass, nowadays, Industries are switching to direct esteri- 
fication instead of transesterification because the former 
route appears to offer some advantages^ over the other route 
and also since- pure fibre-grade TPA is available, Ihe main 
advantage is that the step of manufacturing Ethylene Glycol 
from etl^jrlene oxide can be eliminated since TPA can be directly 
reacted with ethylene oxide , instead of ethylene glycol. It 
should also be noted that there are some disadvantages in using 
Ethylene oxide as the raw material since it is explosive and 
may cause serious safety hazards. In the second stage, i/diich 
is sometimes called the polycondensation stage, BHEI (produced 
from the first stage) is polymerized in a prepolymerization 
reactor upto a of approximately 30 (or upto about 95?^ con- 
version) . In this stage, the overall reaction is operated in 

CJ -f -f 

reaction controlled regime and the viscosity of the 
reaction mass increases approximately to 50 poise. In the 
third stage or in the final stages of polycondensation, the 
polymer from the second stage is further polymerized to a 
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Fig. 1 Different stages involved in the production of PET. 
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of approximately 100 using vacuum and a special agitation 
system. At the end of this stage, the viscosity of the polymer 
melt is of the order of few thousand poise. Due to this high 
viscosity, the diffusion of condensation product becomes 
difficult and hence mass-transfer effects should also be 
accounted, in the overall rate of polsrmerization reaction. 

Under these conditions, it is hecessary to carry out the - 

reaction in special reactors and usually wiped-film reactors^ 
which have facility of mechanical generation of interfacial 
area ace used in the continuous process. In the batch process, 
a single reactor with conventional helical screw or ribbon 

is 

agitators can be used in the final stages of polycondensatipn, 

-jq 

In the fourth stage or the solid state polycondensation stag?, 
the resulting polymer from the final stages of polycondensation 
is solidified and. prepared in the form of small chips. 

Using these chips, the process is carried out at a temperature 
which is in between the glass transition temperature and the 
melting temperature of the polymer. 

Chemistry of DEI formation ; 

In the second stage or the polycondensation stage of 

PEI f ormation, the main reaction talcing place is the poly- 

2 .'': 5 ’ ' 

condensation reaction given by reaction (1) in fable 1,* 

Several side reactions are known to take 'place in addition to 
the above main reaction,^' Ihese are: (1) Ee distribution 
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reaction, (2) Cyclization reaction, (3) Diethylene glycol 
formation, (4) Acid end group formation and (5) 7£rsyl and 
group formation, llie last three reactions are also given in 
Table 1. Due to these side reactions, water (w) and Acetal- 
dehyde (a) are also formed in addition to EG-, Due to a number 
of reactions talcing place simultaneously, the reaction mass 
consists of various functional groups which are also summarized 
in Table 1. The growth of the polymer chains takes place due 
to the main reaction in Table 1 and the quality of the polymer 
formed is strongly dependent on the amount of various side 
products formed. It is therefore, important to control these 
side products,^ 

Mechanism of Catalysis ; 

In the synthesis of PET, the rate constants for the 
poly condensation and degradation reactions are dependent upon 
the aatalyst used in the process. Since the various properties 
of PET are dependent upon the rates of these reactions, the 
selection of a suitable catalyst is important. 

In the transesterification process, the catalysts 
commonly employed are acetates of zinc, manganese, calcium, 

sodium, lead and magnesium. The relative efficiency of these 

2t 

catalysts has been studied ^ and it is found that zinc acetate 
is the most effective. But in Industries, usually a mirfcure 
of these catalysts is used and the catalyst with lower efficiency 
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( SO diata acetate) serves as a ^tracer’ in many cases. In 
the direct esterification process, the proton formed during 
the dissociation of carboxyl group serves as a catalyut. For 
the polyeondensation reaction, the catalyst generally used 
is Antimony trioxide or Antimony acetate. Two tjrpes of catalysts 
are used in the transesterification and polyeondensation 
reactions because these catalysts have their own medium 
dependence. For example metal acetates are very active in 

a high and low lydroxyl content medium but are easily 

25 

poisoned by even small amounts of acid and groups in the 

reaction mass, whereas the activity of antimony trioxide or 
antimony acetate is not affected by acid groups, but it is 
known to increase as the concentration of hydroxyl end groups 
decreases?*^ Ihe mechanism of catalysis is very important 
for a better understanding of the catalyst performance. 

Iransesterification Catalysis : 

The rate of transesterification reaction depends on the 
concentration of a soluble metal alcoholate which is formed 
according to the following reaction: 

M(CHjC0)2 + 2R0H^:^ M(0R)2 + 2CH^C02H CD 

Under certain operating conditions in the reactor, the 
formation of this active alcoholate can be favoured by removing 
the acetic acid, formed in the reaction, fhe rate of trans- 
esterification reaction is reduced if the UW used in the 
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reaction is not pure i.e. if it contains any nonvolatile 
carboxylic aciis such as IPA as an impurity or if there are 
any acid end. groups present in the reaction mass formed as a 
result of side reaction. The mechanism proposed^ for trans- 
esterification reaction employing metal ion catalysis is 
shown in Pig, 2(a). 

Polycondensation Catalysis ; 

in polycondensation catalysis the first step consists of 
the coordination of the metal ion to the ester carbonyl bond 
which enhances the polarity of this bond and thereby facili- 
tates the nucleophilic attack.^^ ihe polarizing action of a 
metal ion is a function of its charge as well as ionic radius 
and the catalytic activity increases if the cationic charge 

increases and the radius decreases. Ihe metal ion can even 

l-'i. f'l ofj 2*7 

build the nucleophilic agent in a complex 'so that its 

attack on the carbonyl carbon atom that has become positive is 

made easier. The polycondensation reaction proceeds by the 

activated hydroxyl group on the carbon of the carbonyl group 

coordinated to the metal as shown in Pig. 2(b), 

Due to the coordination of antimony compounds with 
hydroxyl groups, they cannot react with carbonyl groups of 
esters^ and hence their activity is less in the early stages 
of polycondensation. When the reaction proceeds, the concen- 
tration of the hydroxyl groups decreases and hence the catalyst 





mmrn 

\ / 

O 


/ 


OS 


o 


\ 




sx 


0=0 


o 


o 

(X -cr 


11 


o 


5sr: 


m 

(X 

m, 

f«* 

QC 


a 

o 

CC 


oe I- 

\ / 

o 

% 

’'I. 

DC 


m 

q: 

/ 


2—0 


o' 

! 


IT 


o 

QC 


1 i 


m 

CC 

/ 

□c-o 

mw ' ' 

QC 

\ 

O 2 

\ / V 

0=0 O O 

«* J 

oc CC CC 


c 

_o 

u 

a 

0> 


c 

.2 

#<«# 

o 

u 


■ 4 |m# 

en 

^ m 

c — 

■ O 
u. >* 


O O 
u 

•X? ^ 

0^ o 
01^ ~ 

o 

d. r-i 


Ql 




£ 

.2 S’ 

C •"* 

a ^ 
x: J 2 

O Q. 
<!# C 

"ier %m 

^ Ql 

3 

CM 

m 

Ll 


X 
O 

•>t 

X 
ij 

°V 


? 

O 

(.,J X 

'£ u 






/ % 

(O 

O ' O 
\ / 
u— o u 

X I 

4* 

? 

•*€ 

X 

<i> 

C-5 

I 

u 

fM 

o 

sr 

X 

<M 

u 

tv 

0 
u 

1 

■a 

X 

to 

u 


% 

t 


X 

m 

o 


o 


CM 

X 

o 




\x 

o 

/ 


xo~ 

'«» 

X 

IS 

o 

cS* 

u 

'X 

m 

o 

a 


IS r 

x/ 

o 


0 o 

. \ / , 

1 '.imm I 

1 .^ ^ 

'/ \« 

o o 




u 


f 

to 

u 

X 

o 

cS* 

IM' 

u 

;. , j*4 

o 

"X 


\ 


X 


•o 


/ 


sr 

■u 


Fig. 2(b) MeGhanisin' proposed for polycondensotior^ reaction 
employing metal ion catalysis. 
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activity is improved. By decreasing the melt thickness, 
the efficiency of the removal of volatile component (EG) can 
he improved, thereby improving the availability of the catalyst 
to the polycondensation reaction which increases the rate of 
polycondensation reaction, Tomita"' found that the 
stability constant of dibenzoyl methane complex of each metal 
species is critical in determining the catalytic activity of 
metal compounds in transesterification, polycondensation and 
thermal degradation of PET, Metal compounds having stability 
constants around 12 were found to be more effective as catalysts 
in the polycondensation stage. 

Prom the reaction mechanism proposed, it is obvious that 
during the polymerization there is a molecular weight distri- 
bution of oligomers due to the several side reactions taking 
place in the polycondensation stage of PEI formation in addition 
to the main reaction . Even in the case of first stage reactors 
it has been shown that several side reactions take place along 
with the BHEjO formation. But the amount of these side products 
formed is a very small quantity and it can be safely assumed 
that the feed to the polycondensation stage is BEES! alone even— 
thou^ some recent studies?^*^^ believed that these cannot be 
ne^ected for any true simulation. In order to get the detailed 
molecular weight distribution of oligomers, all the side reactions 
must be taken into consideration and if this is done, the result 
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is a set of lii^ly nonlinear differential equations which 

cannot be easily solved even on a computer. The only way to 

avoid, this computational difficulty is by carrying out the 

analysis of the reactor an terms of functional groups, even 

though it has been shown in the case of PST that one cannot do 

this due to cyclization and redistribution reactions. According 

9—11 

to Komar et al. , if these two side reactions are included, 
the reactor analysis should be carried out in terms of indi- 
vidual polymer molecules instead of functional groups. But 

5 8 

the simulation studies of Ravindranath ' and Mashelkar confirmed 
that the cyclization reaction can be neglected . Also the 
redistribution reaction can be neglected since it affects the 
second and higher moments of functional groups. Hence, in 
our kinetic model, we neglected the cylization and redistri- 
bution reactions. 

In the first stage of PET production, the usual molar 
ratio of DMT and EG in the feed is about 1:2 and the reaction 
occurs mostly in the temperature range of 140-200°C at one 
atmosphere pressure. The use of a catalyst is common so that 
the reaction is completed in a reasonable residence time and" 
it is already mentioned that zinc acetate is the most effective 
catalyst, PED is formed when the product from the first stage 
reactor is heated to 270 - 285®G with continuous evacuation 
to low pressur es^*^ (2-1 mmHg) . The high temperature is 
required in order to maintain the polyester product in molten 
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form, The growth of the polymer molecules is accompanied 
by a substantial increase in the fluid viscosity and extreme 
diffusional resistances are encountered by ethylene glycol (EG) 
and methanol species if the polynner is solidified partially, 
low pressures are therefore necessary to drive the polymeri- 
zation reaction to completion. It is necessary to vaiy the 
temperature and pressure within the reactor as polymerization 
proceeds and the temperature and pressure histories normally 
used in Industry are given by the following equation: 

t = 290 - 50, (a) 

1^3 = 760 exp [ -1.859 0 + 1.776 0 ^ - 1.253 0^] 

(2 

where i: is the temperature in °0, is the total pressure 

in mmHg, and Q is the residence time of the reactor in hours. 
Erom the above equations it is clear that the temperature 
should be continuously increased whereas the pressure decreased 
as the polymerization reaction proceeds. 

Even thou^ PEI is an important fibre from commercial point 
of view, only few optimization studies are available in the 
literature. Based on the kinetic studies of Challa^^ and 
others, Ault and Mellichamp^ optimized the two-stage production 
of PEE. Iheir objective in the first stage was to convert IME 
through ester exchange with EG to BHEI . In the second stage ^ ^ ^ ^ 
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the pressure was reduced to drive off the excess EG-, fhey 
used metal acetate and Antimony trioxide as catalysts in the 
first and secondstages, Ihe reaction temperatures used were of 
the order of 150-200°C. They found that the total reacting 
volume decreased by 25^ during the course of the first stage 
and their studies suggest that the reaction time should he in 
between 25 min and 125 min as equilibrium conditions were 
attained at the latter time. They also recommended a reaction 
time of 70 min. for the second stage. Further studies on the 
optimization of EET reactors include those of Kumac et al.''"^ 

These workers proposed an objective function which would be 
relevant to designing a new plant.^ In the first stage reactor, 
their objective was to maximize the conversion of DMT while 
simultaneously minimizing the amount of various side products 
formed. In the second stage reactor, their objective was to 
bring the of the polymer to a desired value while simulta- 
neously minimizing the amount of various side products formed. 

The proposed objective function included three weighting 
parameters which takes values depending on what is required. 
These weighting parameters were shown to have considerable 
influence on the optimal temperature profiles with some values 
of these weighting parameters, the optimal temperature profiles 
obtained were of the same form as calculated using Eq.(2) . In 
this work, we used the same objective function. 
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In. polycondensation reactions, pressure is an important 
variable besides temperature because the flashing of the 
volatile components, ethylene glycol, water and acetaldehyde, 
depends upon the pressure applied in the reactor. In Industrial 
reactors, both pressure and temperature are independent variables. 
Earlier studies * considered temperature as the control 
variable only and computed optimal temperature profiles using 
various reactor pcessures. We first considered both temperature 
and pressure as the control variables and carried out optimi- 
zation studies for this two-dimensional problem using the 
gradient method. Preliminary computations showed that the 
pressure falls first and settles down on its preset lower 
limit and only after this has happened, the temperature profile . 
begins to be adjusted. These results suggest that the 
polymerization should be carried out at the lowest possible 
reactor pressure. In view of this, in the following we have 
considered reactor temperature as the only control variable. 

While computing, we assumed that the flashing of BGr,W and A 
was occurring at the end of small discrete time intervals with 
polymerization occurring during this time interval. To be 
able to determine the concentrations of volatile components in 
the reaction mass, a vapour-liq.uid equilibrium of Raoult’s law 
type is used. A suitable Hamiltonian has been defined and 
the resulting state and adjoint equations have been solved, 
using fourth— order Itoge-Kutta met^ Using the^^^^^^^ 
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the control vector iteration method for first variation 
tecMique, optimal temperature profile has been determined. 

For the second variation method, the initial guess temperature 
should be very close to the optimal results; otherwise there 
is no convergence. Also the first order technique becomes 
very slow, especially as the optimum is approached. In viev7 
of this, we have adopted the combined method in which we switch 
from the gradient method to the second variation technique 
when the convergence becomes very slow. Due to nonlinearities 
involved in the state and adjoint variable equations and also 
due to memory storage space problems, considerable numerical 
difficulty is encountered in the optimization procedure. Ihe 
storage problem has been overcome by approximating the actual 
temperature profile by 100 piecewise continuous curves and 
obtaining the intermediate values by linear interpolation. The 
results were found not to change when the number of piece-wise 
continuous curves was increased to 200, thus justifying the 
stability of computations. 
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EOHMTMiriOH 

The kinetic model used in this v/ork is given in Table 1 
along with the respective rate and equilibrium constants. 

The rate constants are found out when zinc acetate catalyst 

—4 

of 5*6 X 10 moles/litre concentration was used* Based on 
the reaction mechanism proposed for the poly condensation stage 
of PET formation, the mole balance equations for batch reactors 
can be written and are given in Table 2, In these equations, 
the lower case symbols represent the moles of various functional 
groups and Y the volume of the reaction mass to be defined later. 
The rate constants are functions of temperature and Arrhenius 
relation is used for expressing this dependency. The equilibrium 
constants are found to be independent of temperature. The 
activation energies and the frequency factors of the respective 
rate constants along with the equilibrium constants are 
summarized in Table 3. In the mole balance equations, Z 
represents the total amount of reacted bonds in the reaction 
mass. A separate mole balance equation is not written for Z 
because it is related to the moles of various functional 
groups t hr ou^ stoichiometry. It is given as 

2 - ®g - ®0 - 


( 3 ) 
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glBlE 1 

Reactions InvolTed in the Polycondensation Stage of 

Formation 


Main Reaction 


(1) Eg+Bg 




Z + S 


• Important Side Reactions 

(2) Acetaldehyde Formation 




E, 


g 


^ Ec + A 


(3) Eiethylene G-lyool Formation 

ICt 


E + G 
g 


E^ + E 


E + E 
g g 


(4) Water formation 




E^ + G 


ky/K^. 


^8 


■** 


g 




E^ + w 


Z,4-W 


kg/Kj 


(5) Viny l Gr oup For m ation 

Z 2—^ Eo + 


k^ 


> Z + A 


Eg 

Symbols ; A : GB^GHO 


D : OHGjB 2 OHgO GH2GH2OH 
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OBable 2 


M^le Balance Equations for Yarious Species in Batch 
Reactor (State Variable Equations) 

o 

= ^1 = -al = B3-K4.-E5-2E6+iyE8-RTO) 

. be 

( 2 ) %2 ~ ^2 “ “ 3 ^ ~ + ^ + ^6 ~ ^7 *" ^8 ^ 9 ^ 
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^3 ~ % 

- #= 


o 

(3e 

(4) 

II 

= ■SFT" = *^^9 ” ^10.^ 

(5) 
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(6) 
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^ = y(R^ + Es “ V 


(7) + ®10 - 

^ %) 
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Ej 
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4 
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^6 
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where e^, e^, and e^ are the moles of functional groups 
Eg, Eq, and E^ at time t and that of Eg at t=^. 

The mole balance equations in fable 2 represents the 
equations for the state variables X^, X 2 , X^,X^.-*X^ and are 
written in a compact form as 

0 dX 

S = -ffE = f (X, I) (4) 

where X ■= |^ X^, X2, X^, ] (a) 

£ “■ L » ^ 2 .* ^3* ^4*^ ^ (h) 

(5) 

Since the feed to the polycondensation stage is assumed to 
be BEET alone, the ’initial* conditioips for the state variables 
can be given as 


x^ (t=o) = 1.0 

(a) 

3^(t=o) = 0.0, i = 2,3,4^, 

a) 


(6) 


fhe molw balance equations in fable 2 cannot be solved 
without knowing the quantities Qj;q, % and which are 
present in the mole balance equations for the volatile components, 
equations (5), (6), and (7) in fable 2. fhese quantities stand 
for the amounts of ethylene ^ycol, water, and acetaldehyde 
flashed from the reaction mass due to the boiling phenomena 
occurring at the temperature and pressure applied on it. In 
order to de'fesrmine these unknown quantities, a sound theoretical 
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mpdel should be developed which relates the amounts flashed 
to the conditions existing in the reactor i.e., to the 
temperature, pressure, and composition of the liquid boiling. 
Here, we proposed a simplified model for the flashing process 
in which it is assumed that the flashing is occurring at the 
end of discrete times lAhereas the polsnnerization occurs 
without flashing during this time interval. This means that 
QjXj? should be replaced by QgQ_ 6 (-t- - nA t), 

(t - n At) and 6 (t - n A -fe) where ’n’ is the number 
of discrete intervals in time, if. At is the width of the 
time interval and 6 is the dirac delta function. Using the 
above model, we carried out the computations hoping that the 
results obtained would simulate the flashing process. 

Since the polymer molecules are assumed to be involatile, 
the total pressure in the reactor should be equal tothe sum of 
the partial pressures exerted by EG, ¥, and A. Assuming 
that the vapour leaving the boiling liquid is in thermodynamic 
equilibrium With the residual licpi id and assuming Eaoult*s 
law is applicable , we can write , 

Ej = Xj + pX + \ (7) 

where Pj is the total pressure in the reactor, P^, and 

are the pure-component vapour pressures of EG, ¥, and. A 
respectively at the temperature existing in the reactor and X^, 
and X^ are the corresponding mole fractions of the coB^onents 
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in the liquid phase Usually the temperature of the reaction 
mass is greater than the critical temperature of acetaldehyde and 
therefore it cannot exist in the reaction mass. 2^, 
in the reaction mass are obtained from the following equations 


Xg = — i — 

^ (Sg + + d + e^) + g + w 






w 

4- d + e^) “f g 4- w ar 


377 


■& 

d ^ ^ -K W +' S' 


(b) 


J^X^J 

( 8 ) 


where the first term in the denominator in the above equations 
represents the moles of polymer, Mp, Since the vapour 
pressures of EG-,¥, and A ^an be obtained from any standard book 
on theimodsmamic Xq, and 3^ in the reaction mass can 
easily be computed. 

Ihe volume of the reaction mass, V, is computed using 
the relation 

Y = Vp Mp + 'V’g,Xg_ + X-j^j. (9) 

where Vp, Yg., and Y^ are the molar volumes of the polymer, 
ethylene glycol, and water respectively. Yp is assumed to be 
the molar volume of M. Yp, Yq and are temperature dependent 
and the equations expressing this dependency are given later. 

Since our main goal is to determine the optimal temperature 
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profile in the reactor, we should select some method of 
optimization. In this work, we applied the variational method 
of optimization to PET reactor. We assumed the following 
functional relationship for the volatile components, EG,W gad' 

A: 

g = g(T, eg, e^, d, e^; E vu (a) 

w = w(T, eg, e^, d, ' (b) 

ar^ a (I, e^r e-) 

( 10 ) 

and defined the state variable vector as in Eq.5(a). 

Objective Punction 

As mentioned earlier, the quality of the polymer 
produced is largely dependent on the various side products 
formed. Since the reaction time and the pressure inside the 
reactor are specified, the only variable that can be controlled 
is temperature, The temperature should be controlled in such 
a way that the nimber average molecular weight of the polymer 
produced reaches a desired value, in the shortest time 

possible and at the same time minimizing the amount of side 
products formed. It is also found that an excess of diethylene 
glycol, e^, in the final product is undesirable because it 
lowers the melting point of the fibre, but to improve the dyeing 
ability of the polymer, some small amount of if (say E*) is 
needed. Bearing all these objectives in mind, the objective 
function is written as 
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I {[E(-b)} 


a 

o^.(e -D*) ].. + / \ 

1 3 } ^ 4 ^^, 


®g 


+ “3(6^ + e2)}at (11) 

where 1(1;) is the temperature as a function of time in a 
batch reactor to be determined optimally so that I is minimized, 
t^ is the reaction time, is the desired value of the 

number average molecular weight of the polymer, D* is the 
desired amount of diethylene glycol in the polymer and ^2 

and are the weighting parameters representing the relative 
importance of diethylene glycol in the product, number average 
molecular weight of the product and the amount of side products 
formed. It is found that the weighting parameters have a 
considerable influence on the optimal temperature profile to 
be determined. The objective function proposed allows more 
flexibility in reactor design and forces the product from a 
batch reactor to be as close to the desired value, ^^^d' 
possible. This is because at the design stage, a few percent 
deviation in the value of can be tolerated without 
causing serious damage to the physical properties of the 
polymer. 
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IffigHOD OF SOniTM 

With "the definition of the state variable vector in 
Eq..5(a) and. the objective function in Eq.(1l), the optimal 
temperature profiles in a batch reactor can be easily obtained. 
One defines the Hamiltonian, H. as 


e 


'^o « 


\ (®c ®v ^ ^ \ 

i=1 


where are the adjoint variables given by 


( 12 ) 


-d 

dt 


dH 


= 1, 2,5,4 


(13) 


with the final conditions as 


=“fTr [«i(ep--^*) ] 1 ^ i = 1f2,5,4 (14) 

1 t-tf 

Ihe details d equations for the adjoint variables along with 
the final conditions are given in Table 4« 

Having solved the state variable equations. in Table 2, and 
ibe adjoint variable equations, Bq.(13), the optimal temperature 
profile is obtained by using the necessary condition. 


^ = 0 = 


df^ (X^, X^, x^, X^, T) 


(15) 
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^ % 

The detailed equations for a^^e given in Table 5» 

In obtaining the optimal temperature profiles, the 

u. 

temperature is constrained to irie between two limiting 
values 

400°K < T < 600°K (16) 

and it can be shown that the optimality relations described 
earlier can be used for this constrained problem. 

In order to obtain the optimal temperature profiles, 
use can be made of several iterative procedures existing in the 
literature. In all those iterative procedures, at least one 
set of equations from the governing equations for the optimi- 
zations feftjm the gov^ening e-quetiens for tire^optimizarfeioH: 
problem should be assumed, for example, this can be the unknown 
initial conditions on the adjoint variables or the entire set 
of control variables. With the assumed values of any of these 
variables, the equations governing the problem arc solved and 
later, the original assumption is corrected by using some 
technique. The new values calculated yield a result closeTto 
the optimal values than the assumed values. This procedure is 
repeated assuming convergence occurs. In this way at least a 
local optimum can be found out. 

' lynamic programming can be used as a method for the 
calculation of optimal temperature profile. But this method 
has a drawback in that it needs a eonsiderable large memory 
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storage space for sol-ving system equations in which there 
are only three or more state variables. Despite this drawback, 
it can be used for systems in which there are constraints on 
control and state variables. Inother method called the 
"gradient method" or "control-vector iteration procedure" 
as suggested by Eay and Szekcly^^ can be used. In this method 
the system equations and the boundary conditions are fixed 
and an iterative procedure developed which converges to the 
desired performance index, 

¥o have chosen the control vector iteration procedure 
for computing optimum temperature profiles. Moreover, the 
algorithm for this method is simple to program for a computer 
But the disadvantage with this method is that it becomes very 
slow especially when the optimum is approached and does not 
converge to the true optimum. At this stage the computations 
should be stopped and some other method must be used for 

using the convergence or for finding the true optimum. As 

57 

suggested by Dapidus and luus , we used the second, variation 
method for achieving the true optimum. But this method is found 
to lead to a lot of instability problems if an incorrect 
initial guess is used. In order to overcome this problem, we 
used a combined method for solving our optimization problem. 

In this, gradient method is used for the first few iterations 
or till the point where the change in the value of the 
obyecti’^^e function from iteration to iteration is small anil 



at this stage, one switches over to the second variation 
method until the optimal profile is obtained. We found that the 
combined method works better than each method applied 
separately. These two methods are described separately below: 

Control Vector Iteration Pcocedure : 

In this method (sometimes called the gradient method), 
one assumes a temperature profile T^Ct) . With this the state 
variable equations (Table 2) are integrated in the forward, 
direction (from t=0 to t^) , storing the values at small time 
intervals. The objective function is computed using these 
stored, values of state variables. The adjoint variable equations 
(Table 4) are integrated in the backward direction from t=t^ 
to O) and finally the assumed temperature profile is corrected 
using the equation 

TneW(-t;) = T°^‘^(t) - (17) 

where e is a parameter which determines the magnitude of the 

Qorrection to be employed in the old temperature profile and 

is assumed to be independent of time t. Assuming a value of e, 

the temperature profile is corrected and the above procedure 

is repeated till the necessary condition given in Eq.(l5) is 

' B H 

satisfied. It is important to note that the value of in 
Bq,(l5| in the first iteration will not necessarily be zero 
since the assumed temperature profile is not optimal. This 
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method is extremely efficient and the optimal temperature 
profile is usually obtained after two or three iterations* 

In the literature, there exist several techniques for 
finding the value of e from iteration to iteration. In this 
work, we selected that technique in which an interactive 
computer terminal is used. In this technique, various values 
of £ are assumed and the new temperature profiles are obtained 
using Eq.(l7) with the temperature profiles obtained, the 
integration of state variable equations is carried out and 
the corresponding objective functions are c^culated. The 
various values of e and the corresponding objective functions 
are examined on the terminal and the approximate value of ^ 
corresponding to the lowest value of the objective function is 
selected as the optimum s for the next iteration of calculations. 
With the value of optimum e , the old temperature profile is 
corrected and the entire procedure is repeated till the value 
of the objective function does not change significantly from 
iteration to iteration. The technique suggested by lenr^®, 
does not require an interactive computer terminal, but there 
are several situations when this reads to unstable results. 

We selected the one requiring since it yields the optimum e 
in a minimum number of trials. 
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Integration of State Variable and Ad .joint Variable Equations ; 

The integration of these equations is carried out using 
fourth-order Rung e-T?:utt a method. As is evident, the integration 
of the state variable equations (Table 2) is not possible 
since the quantities QgQ.j and are unknovni. However, 
to overcome this problem, we. used the method proposed earlier in 
this thesis. In this method, for a given time interval, the 
integration of the mole balance equations for the volatile 
components is carried out without the terms and 

and then at the end of this time interval a flash subroutine 
is used to find the concentrations of ethylene gLycol and 
water such that the Eq,(7) remains satisfied. In the flash 
subroutine, with the values of the state variables obtained after 
integration time interval At, the total pressure Pq, is computed 
using Eq.(7) and if the computed value is found to be greater 
than the actual pressure applied in the reactor, then it is 
obvious that some amount of flashing has occurred. The amount 
flashed can be found out using equilibrium flash calculations. 

If the computed pressure is found to be less than the actual 
pressure, then there is no flashing. This way, the state 
variable equations are integrated for the entire time of 
polymerization. The total time of polymerization has been 
assumed to be two hours and this time has been divided into 
2000 equal intervals in order to avoid excessive memory storage 
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space. But when the o-omputations are carried out with a At 

—3 

of 1x10 hr, a numerical instability was encountered., because 
of the fast changes in the concentration of EG-. It is also 
found to become negative with At =10"’^ hr, when the vigorous 
flashing of ethylene glycol starts at lower pressures. We 
found that by decreasing At to 2x10“^ hr, numerically stable 
results could be obtained. But vdaen At is reduced to this 
value , there is a problem of memory storage . However, we have 
overcome this problem, by storing the values of state variables 
at every time interval of 10“^ hr or at 2001 points of time. 

The integration of the adjoint variable equations (Table 4) 
is carried out with At of 2x10“^ hr, usiag linear interpolation 
for the intermediate values of state variables. The adjoint 
variables are also stored at 2001 points of time. 

Eor all the above integrations to proceed, some initial 
gueffi should be used for the temperature profile. Here we 
assumed that the temperature profile is isothermal and approxi- 
mated the actual profile by 100 piece-wise continuous curves. 

In the integrations, the intermediate values of temperature 
are simply found by linear interpolation. 

A listing of the computer program for the control vector 
iteration method is provided in Appendix. 
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Second Yariation Method: 


As mentioned earlier, this method I'eads to a lot of 
numerical instability problems if an incorrect initial guess 
for the temperature profile is used. In order to overcome these 
problems, control vector iteration method is used for the first 
few iterations and with the near optimal temperature profile 
obtained using this method, the present method isstarted. 

In this method, with the near optimal temperature profile 
obtained from the control vector iteration method, the state- 
variable equations (Table 2) are integrated in the forward 
direction (from t=0 to t^) storing, the values of these at every 
time interval of 1 z 10”^ hr. The objective function given 
by Eq.(l1) is computed using the stored values of state variables 
and the adjoint variable equations (Table 4) are integrated in the 
backward direction (from t = t^ to O) . 3n addition to these, 
the following equations are integrated in the backward direction 
(from t = t^ to O), 
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n J* 

The detailed equationa for the various terms, ^ ■ , 
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equations (18)-(29) are given in Tables 6, 7, and 8 respectively, 

afi 

The equations for — rffi— are already given in Table 5. 


In Eq,(18) 
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The detailed equations for the elements involved in the above 
matrix are included in Table 4 . 

Equation (18) is also known as the first— order Siccati equation. 
In this P is a symmetric matrix of the order (4x4) as given in 
Eq.19(a) , It is also related to the adjoint variables by the 
equation 

P = (4^ (31) 

The final conditions on P are: 
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.2 

P(t^) =_5_2. [a^ (e^-D*)] (32) 

6 X 

In Eq.(27), £ is a vector function of the order (4x1) given 
by Eq, 28(a), The final conditions on £ are 

£ (t^) =0 (33) 

During the backward integration of the equations (18) and (27)» 

9 f ^ . TT _2 tt 

the values of P, £, R, — " - -- w - ' are stored in ■ 

- - O-L -L 

the computer memory. Using these stored values, the pert'urbation 
equations for the state variables given by 
d( 6 X) 9 f^ 9f^ 

5 = (?ir) <5X + ^54) 


where 


iX = [dX^, SXg, 6X^, 6X^ ]■ 


d 5X. d 6X^ dSX^ ddX. T 

* - ^ ^ - J- 1 f'hV 

» dt * dt » dt J 


(35) 


are integrated in the forward direction (from t=^^ to t^) with the 
initial conditions as 

6X (t=0) =0 (36) 

The above initial conditions are obtained noting that the initial 
eonditions on the state variables given by Eq.(6) are constant. 


The values of 6X are also stored in the computer memory 
. — • 3 f , , 2 

With the stored values of P, £, R, » (-|^) » ( — -"g )y and 

the ten^je nature profile is corrected using the equation 




42 


0 -L 9j 

“'[(— 2 )"**'^ (37) 

where 'j* stainhs for the iteration nunher and e is a step 
size parameter determining the magnitude of the correction 
employed in each iteration, This parameter is introduced hy 
Merriam to prevent overstepping in the temperature correction 
and also to maintain the validity of the approximations used 
in deriving the Eq,(37)« According to Merriam, the limits set for 
e is 0<e<_1, The value of e is found out using an inter- 
active computer terminal. Eor the first iteration, the value 
of e is taken to be one and with the temperature profile 
obtained from Eq.(37)j the state variable equations are integrated ; 
and simultaneously the objective function is computed using 
Eq,(11), If this objective function is found to be greater than [ 
the previous one, the value of e is halved and the iteration 
is continued, This procedure is repeated' till the objective 
function decreases. The value of e corresponding to this 
objective function is taken to be the optimum e and the 
temperature profile is corrected. This completes one iteration 
of computation. Ihe entire procedure is repeated till ^ he 
objective function does not change significiantly from iteration 
to iteration. It is found that with a value of e =1.0, faster 
Gonvergenee is attained i.e. the lowest value of objective 
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function is obtained, and if e > 1.0, there is no convergence. 

She optimal temperature profile is obtained using the 
necessary condition given by Eq.(15) and the sufficient conditions 
given by the equations 

c? H 

— rr~ > 0 for 0 < t < t^? (a) 

]P(t) finite for 0 <t 1. t^ (b) 

(38) 

According to Bryson and the above sufficient conditions 

are knovna as (i) the convexity condition (or strengthened 
legendre-Clebsch condition), and (ii) the condition that no 
conjugate points axist on the path (or Jacobi condition). Ihese 
two sufficient conditions determine whether any neighbouring 
extremal paths exist. If ary of these sufficient conditions is 
not satisfied , then the second variation method will diverge 
which in turn iirplies that the neighbouring extremal paths does 
not exist. If P(t)^ oo at t=t*, where 0<t* < t^, the 
integration of Eq.(2(!r) should be stopped because there is a 
conjugate path for t< t’ , This behaviour of P, we encountered 
in our problem ard as suggested by Bryson and Ho, we stopped the 
integration at this point and assumed that the temperature 
profile for t£ t’ does not change or it is the optimal profile 
that can be obtained under the existing conditions. By using 
this cohcept, we have obtained the optimal temperature profiles. 
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The integration of this equation is carried out using a 
fourth-order Sunge-Kutta method. iTrom the computational point 
of view, this is a very important equation since almost all the 
variables involved in Eq.(37) are obtained by solving this 
equation. The integration poses both a computational as well 
as a memory storage problem becai:©e of the following reasons: 

As in the integration of the state and adjoint variables, if 
a step size of 2x10~'^ hr is used, there is a problem of memory 
storage since P is a (43:4) matrix and it is difficult to store 

the 16 elements at 5000 points, of time. If the step size is 

—3 

increased to 1x10 hr or more, then there is a problem of 
numerical instability and the integration using fourth-order 


Eunge Kutta method could not be carried out. In erder to overcome 

—4 

both of the above problems, we used a At of 2x10 hr in the 


integrations, but stored the values of P at every time interval 

—2 

of 1x10 hr i.e, at iCfDpoints of time. The values of S, 



i—o): 


dfi 

— are also stored at 101 points of time. 


The stored values of state and adjoint variables at 2001 points 


of time are used in the above integrations, and the intermediate 


values of these are obtained by linear interpolation. In this 
way we would get stable solutions. One important point to be 
mentioned is that,' in the integration of this eqn., whaxever 
P oo for any t* , where 0 < tV < t^, we stopped the integration 
at this point and assumed a Talue of zero for all the elements 
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of P at lespar times i.e. for 0 <, t <_ t*. 

Integration of the q. Equations : 

The integration of these equations are carried out using 
a fourth-order Runge-Kutta method. Unlike the Riccati equation, 
this equation poses only a memory storage problem since £ is 

a (4x1) vector function of time. In order to overcome this 

—5 

problem, in the integrations, At of 1x10 hr is used and the 

—2 

4 elements of q are stored at every time interval of 1x10 hr. 
V7e could get stable solutions using the above step size and 
also the main advantage of using this step size is that the 
values of state and adjoint variables, available in the memory 
at every time interval of 1x10 hr can be directly used in the 
integrations without any linear interpolation. 

Integration of the Perturbation Equations of State Variables : 

These equations result from the linearization (Taylor- 
series expansions to first order) of the non-linear state 
variable equations and are given in Eq.(54) with the initial 
conditions given in Eq.(36). The integration of these equations 
is difficult since both 6X and 6T are unknown quantities. In 
order to find 5X, we should know 6T before-hand but since 
'^ur main problem is to correct the ten^jerature profile using 
6T is also •unknown and hence calculation of 6X is 
not possible. Normally, to find 6X, we have to solve a 
two-p©int boundary value problem formed by the equations 



46 


6| = A(t) 5X-B(t) 6^ 

6X = -C(t) 6X - A®(t) 6X 


■where : 

A(t) 

B(t) 

G(t) 


df 

df 

(— ) C 


6f 2jr ^ ^2rj 

) r_l_£ 'v V s h: \ 


ai 


as 


aT 

- (-1^) (4'^ 


sfs 

63 ) 2 ^ 


dl 

df 


T 


■ariT^ 


33} 


33} dX 


) 


(a) 


(b) 

(c) 


(39) 

(40) 


(41) 


3}he solution of the above two-point boundary value problem is 
time-consuming as we have to first eliminate either 6X or 5\ 
from ^ne of the above equations using some simplified assumption. 
As suggested by Biyson and Ho, we used the backward sweep 
method in which the solution for 6 A is assumed to be of the 
f crm 

5 A = B 6X (42) 

where P is obtained by solving the Siccati equation. Since 
the value of 6X at t=0 is known, the value of 6 A at t=0 can 
be found out and then the equations (39) and (40) can be 
integrated forward as an initial-value problem. Since this 
is a very tedious process, another way of doing this is to 
simply substitute Eq.(42) into Eq.(39) and then integrate 
B9*(;39) forward as an initial— value problem. But we found that 
this procedure worked only when there are perturbations in the 
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initial conditions of the state variables. Hence the method 
for solving 6X we finally adopted is this, in which the 
expression for 61 (= given by Eq.(37) is simply 

substituted in Eq.(39) and then integrated this equation 
forward as an initial value problem. Ihe values of 6 X so 
obtained are substituted back into Eq.(37) and the actual 
correction in the temperature profile is obtained. Ihis kind 
of method worked for our problem since the values of 6X obtained 
were of the order of 10"* to 10"* , which are consistent for 
assuring neighbouring extremal paths . 

A listing of the computer program for the second variation 
method is provided in Appendix. 



TABLE 6 


DETAILED EQUATIOI'iS EOR THE TERI^S 


a'l. 

az^ 


^"rr 

3X^ 


J 

ax^ 


2t 

a J 


ae. 


■ 2 _ 2 2 
■ 3^ J 3 I 3 I 


3^ J 


asg asgae^ ae^ 3d ae^ a 


2 2 

3 J. 3 J 


3^ J 


3^ J 


ae^SCg 3e^2 ai^d ae^ 3e^ 


-2t 
3 J 


a^J 


3^J 


33 ■ 3e^ 3d3e" 


33 


T 


a2j 


ad 3e, 


2t 

3 J, 


a®-y a®g 


d^J 


9^ £ 


3 e^ae^ ee^ 3 d 


3^ J 


ae^2 


2 So.. 3 e„_ 

2 g*_ f go 


2M. 




eg" 


nd 

T“ 


I 

ae^ 


2 T 

a jL 


ri B. 


T = 2“3 


a"f. 


1 • / \ /9 

■Tj - This is a (4'2:4) matrix of the same form as (.- 


The remaiiiing elements are zero 

.2 

3 

3X" 9X^ 

with ^ replaced every^vhere by f^. The detailed 
equatiore for the elements involved in this matrix 
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fie 2 
g 


-( 


4 k , 

ir ' 

2 kr 


“ T " 


( 


8 k , 4 kp - 

TTj ^ 36 ^ 


s 


ii iCrj 2 kp 


a^f 


1 


89 (, 8 eg 


4 k ^ 2 kp . 

_ ( ^ 4 . j iP ) ^g 

T K^ + “nr ^ ne^. 


( 


4 k , 2 li :„ ^ 

■ ,/ JL ) 

“ T -^ oeg 


ky ^8 \ 3w ^8 / ^ 1 fiw 

^ vk ;;; + tt 5 ) ie 


fiw ■ \ 

■0 ^ % jeg > 


fi ^ f . 

ad fieg - 0.0 


24 * 

9 f^ 

4 kp, 2 

, z' 5 , 

^7 ) 9 g 

■ 4 k^ 

0 % 3®0 

' V ^ 

“T" fie^ 

" 


^8 /fiw 
" de^ 

+ iW_N 


_ 

2 

9 f^ 

9^f^ 

fi^f^ 

9egfia 

geT 9d 

0 

fid^ 

9e^ fid 

_ 

- 96^2 

. ^ 

■ " wr 

2 kg 

“ -w^ ' 

9W 

9'«v 


fie. 


= 0,0 


The remaining elements are obtained by noting that 
■1 


2 4 > 

3 £. 




9 ^ f . 

file . "fiX where ij3 ~ ”5>2j3>4 
31 


3 ^ £ 


2 . 

- ■■ /n - This is a (4x4) matrix of the same form as ( — — ^ ; 

3X dx 


with 
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i replaced everywhere by The detailed equations for 

the elements in'vol'ved in this ma,tri 3 c are; 


a ^ f 2 

■ 

.2 

9®c 9®g 


4 

-T” 

-M. + 

ae^ 

g 

2^6 

■~T" 

2 

+ T 

kf^ 

aw 

2kQ- 

aw 

00 

g 

2 k^ 

Y 

-M. _ 
9®c 

2 ky 

“T" 

a eg 

+ ( 

■®4 ■ 

■'S ) 

aw 

S®c 


8 . 1 


(1 + 

^ ^ Kp - de „^ 

5 g 


3 ^ f , 


= 0.0 


g 


2 

a fg 

i !^2 

3 ®^ “ 


2 k 


5 . ^ ^8) 

Be ^ '^ TTET ^ “ w-J 


4 

- ir ~ 


V ^ 

2 k, 


V 


de 


g 


ae . 


^8 


aw 

3e_ 


a^fp 

= 0.0 

-2 

5 f Q 


6% 


2 krr ko 

-Y-^ ^ ^ a _ w ^ 

V a e„ ^ 5e^ a ^ 


2 

d ^ f . 


'2 a^fp ■ a^fp a^ fp 

a^g 3 'i "■ a a - = PTaT = 


a ^ f , 


de ^2 


^ ^ aw 


Eae remaining elements are obtained by noting that 
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2 f 
d 2 

“aTTlTT 


df 


2 where i,j = 1,2, 3,4 




"^3 

^"„2 ~ This is a (4x4) mairix of ihe sanie form as (- — 2) with 

dX^ 

% replaced everyv.^here by f^. The deiailod equations for the 
elements involved in this matrix are; 

4 ktr _ „ 2 k. 


^ f 3 


aeg2 

=^. 6 M . 

+ _ 

9 g 

b 

Y 

1 

K\ 

e^f^ 2 k^ 


^0 

aepVe^ T- 


d f 


6d 6e 


= 0.0 


'g 


d f . 


2-f 

5 ^3 


2 k 


j_s 
5®v 


9®v ^®g SOg ae^ 

The remaining elements are zero 


Y - This is a (4x4) matrix of the same fomr as ( — -^) with 


J replaced everjn/rhere by f ^ . The detailed equations f»r the 
elements involved in this matrix are: 

k. 


d^f, 


.2^ 

a -i- 


ae^ ae 


4 


ae^r ae. 


3 

“7 


The remaining elements are zero. 
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gABIS 7 

DETAZLSI) SQUAD lOlTS SOS THE TEBJ^SS 


a^J 


dT ax 


io & x 


9^i. 

‘aTTT 


a^J 


a^J 


a^j 


0 J 


^ aTae^" orTd’ TFre“J 

^ G V 

All the Glements in the above matrix are zero. 


iT 




8% 


•j 0 tJ 

- This is a (4x1 ) vector of the same form as (vsr^) 


^ -4_w _l.k> ^ UX4.^ OCJJLi-LO J-.WJ_Xii CL ^ 

with J replaced evcr 3 /where by f -j . The detailed equations for the 


elements involved in this vector are: 

2. 

^ / s ■ 4 a^ a(k„/V) 

- (-4 6g - e^) — ^ [ (_^) -H + g 


■1 

WT^ 




^ ^®go"'®g-®c“®v^ r/^ 3 \ 3^ g ag ‘^^^3/’^) 

+ ^ — „ L (^) a¥Teg+ ^ -W~ 


d (IC 4 ) 


k 


dT 


_ 2e r (^) 9-1^ 

L V'' a T ao 


a T ae. 


+ 


dCk^/Y) 


ae. 


k 

- 2 [ (^) -^ + g 


d(k./Y)_ 

- 4 e 


dT 

d(kg/Y) 

““dT 


-] 


+ 2 e^[(^) rj A fg-- + 


ag 

a'l asr- se. 


aCkyA) ^7, 02, 


+ 


a ~w- 
a eg 

dCkg/Y) 
®c dT" 




dr~"] “ l| c ("T^ 

dCl^^/Y) 


w 


k 




'8\ aw 
Y^ a^ 


'[ (“tt) + 'W 


dT 

dCkg/Y) 


J 


dT 


1 


2. 

8 f 




- if + 8 + 


4 (e — e *-6 — e ) 
^ go g c 

“ — | r"— — — • — - 



53 




d(k^/Y) k„ 2 d(k„/Y) 

+ -gj ]+ 2e^ [ (-^) J 


+ 2 [ (-^) + 
d(k^/Y) 




a V 

ae^ 


-g^- ] - G 


~W 

aCkg/Y) 
gj— 


k 


- ("T^ -fl + w 


a(V'^) 


w 


+ 


(®gO-«g-®o-®Y) 


-JiL-TZl r (; 8 n _2_w 

TZ L ^7“^ aFae ^ ^ 


a(kg/Y) 


W 


■] 


■] 


3^fl 

wn = 0,0 

9 ^1 4 _ /^5 V a CT a(k^/Y) 

“ “ 13 ^ (-4^ -^ + S — “] 


aT 


■'v 


+ 


I5."g0-V°c-^T) r _a!s + (^) 

—■ 7 - “ aiae^^ ^ dT -1 


kp- 2 

« 2e r (~) -■a.r.g, 

g L ^ Y^ 6 5) a e. 


+ 


V 


d(k^/Y) 

C a g... ) _,§„-] 

^ a e ^ d T -I 


2 ^ „ dCk^/Y) 

+ 2e [ — + ("1^ ) ~”gm“~“] 

® '' 3 T a ^ ®v 


e_ trj 2 
(\ aw 


-ecH) 


- Y 8T ae^ 


+ Tnr- 5J i 




^ [ (^) -f-? + W 


aCicg/T)^ , <®so-V®o-®t) 
31 Ko. 


[(4) 


.2 „ ,, d(k„/Y)' d(k,/Y) 

_ aw r h W s _.2l__J c — -^3^ 1 

3T ae ^ae^-'^ ® ^ ~ ^ 


"g d!E 
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" 5)his is a (4x1) vector of the same fom as ( g-y q y ) ^4th 
J replaced everyi-/here by ± 2 * detailed equations for the 


elements involved in this vector are: 


<3(k/) 


dl de. 


+ 2 e r(~) 

+ aT 5' 




3 / ^ / 

1 

m J 




d(k^/7) 


2 


kv .2^ _ dCk^/Y) 

“7^ “"dl ^ 

a S 


+ iq- aT ae + w: — 3T ■' 


+ 1_ ff-l) -M + 
~ “aT ^ 


d(ky/Y);:, d(kg/Y) 

^ di ' ^ ” ®c _ p 

a(V^4 , (°gW%-^Tll8 


J - e. 


(-n) 


+ -Lw 


dCkg/Y), 


^ d(kg) 


^ "^2 „ „ 2 q r (!i) 3^ g + 

did 6^ Y^ al ae^ . ae_ dl 

C/ (JO 

ICy 2 aCky/T) 

-2e f f—w) 3-^ + — -4fl — J 

al ae^ 3e^ . 


■ k^ 
-2[ (^) 


d(ky/Y) _ 


"U 3 ^^ 

J+ ^ L(“ir^ 


d(kr^/Y) - 
+ ^-sr-'— ] -e 


d(ko/v) 


> 2^ r(^) 


3\ a w 

Y^ 3 1 


a(kg/T); (O -e -Q^-op kg 

— frp — ]+ L'-t’ 
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-) aCk^) 

^ -w 


aT ad “ * 

9% <3(k /V) 

dT ao^ g^W^aiae^^ ae^ ai’ J 


itry 2 d(kr7/V) 

■2 So [ (- 7 ) 3T%g- + —m — ^ 


a d ac^ 9 Gy 


®grr^7x a^w 


+ ^ [ (-4) 


aw 


ad ae^ ao^ 


ai J 


1 r / 8>, aw 

" ^ t C-7) -^ + W 


a(k„/v) 


2h:- 

9 ^ 
id ax 


' 5 

ad oe. 


.'^gu-V°o-Vs r ,'^S^ {,2^ , 0W ■ 1 ^(kg) 

^ ^ ) !_ ^ "“■3d — 3“ 7 *~d[!r 

5 


.2 J 

dhis is a (4x1) vector of the same form as ( qY "a ~ x 
with J replaced everywhere by f^. dhe detailed 
equations for the elements involved in this vector 


are : 


2 r (-^) 9._1£ . LjS 

eg ad 3G„ ^ ae^ 


d(kc/v) 


dd 4 


dCk^/Y) 


+ 2 [ (-^) +ng 4- 2e 


adr^/T) 


flg aCkjA) 

ad ae^ g Y^ ad ae^ ^ ae^ ad J 


2 

a f. 


= 0.0 
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dTTe^, 


A 2 

2 Sg [(-1^) a 


4!^ 

d 2 d e^ a e^ ai 


A 

dTlT 


T 

This is a (43:1 ) vector of the same form as ("^jr-^) 


with J repi.aced, everywhere by f . . The detail ed. 
equations for the elements involved in this vector 
are: 


2 

Wae. 


2 f ■ 

a ^4 


1 


aCkg) 

■~a!r 


dT a e. 


1 


- e 


V -“df-— 


ClCkg) 

dF” 


a^f 


eff 3d 

2^ 


4 _ 


= 0.0 


a "-4 

1 

- ^ 

""W“ 

- "g 

3 ^ D 

, 2 

2 



3 g 

_ 3 

g„.= _ 

aFael ^ 

ai 3 q^ 

ai 

6 d “ 


dCk^V) 

'w 


dlae. 


2 

a "w 

a'i' ae. 


2 2 2 

3^ — d — d 

dFle. a'i' ad 


■T 


2{P»-Pl,) 


■T 


31 ae 2®;^-%) 


i!? 

dT 


dP* 
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2AB321 8 


lETAILED EQUATIOFS POR THE TEEMS 




- 2 e 2 . 1±So;;!£V:!i 1 r & ^ 


— &-■ 1? y_- T s. 

.m2 


d^(k,/v) 

+ g 2“ — + 2 

dT^ 


d(k,/v) 

y ^ a gi 

-~m — af J 


d^(k^) 


k.^ 2^ a^(km/v) d(k./v) 


d^Ck./Y) 


+ 2-1^^ -^1 
^ dT dT-i 


kr, .2, 


d^(km/v) 






d(k/Y) ■ 

^ dT ^ ®c®g 


d (kg/Y) (eg^-eg-e^-e^ 


^2 + ^ 
V 0 


d^(kg/Y) d(kg/Y) 

ai2 + 2 55- 


aw . 
“dT J 


•e-rr^n. 
V g 


d^Ck^/Y) 


f. 


d^(k 


^ + 2 e f(^5) + 

2 + g*-^ Y^ 31 ^ dT ST ^ 


dCk^/Y) 


d^(km/Y) _ 


J + ®p-2 


d‘=^(kg/Y) 


d^Ckm/Y) 

" -^ + S — V~ 

° dT'^ dT'^ 
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, ®g r (^ 7 ) 3^w 
4 33) 


a^Ck-^/v) a(k^/Y) 




ill 


^ _ a^kg/V) (e -e-o^-e^) kg 2„ 

■-^o ^ [(^) ^ + 

„a2(kg/Y) aCkg/V) ^ (e -e -e„-e^) a^Ckg) 
—^+2 ^ 


% o /’"Sv a^K a (kg/Y) 

— r = 2eg [(-4) 2-Sj + g 2 + 2 

83 ) 6 i. V ^^2; 

+ e 2 


aCkj/Y) 

a3) 33) 


_ (fgo-®g-®o-®Y) (*^9) 


a^(k,/Y) 


- e e 

T g 


o 

a p° 

G 


2 

I"® 

w 


3 


- ^ ^ 


a^ p° 

-g ^ - w 


- Pg) 

a^ p° ap° 

r W o w 3 W 

~ ^ “ 3 !r "TIF 




a^s 


|_ = iaOii] +[2229 . 1^] 


"G L 5 


2 

a p° 

w 


po r 3515.706 r 1757.853 ' ] ^ 

^ , (P-33. 274)^ (3)-33.274)^ 


a^(k^/V) 


( 3 !ki* . E* m 2 <^17 -% 

i -I - (21 + ||) 5 



59 


1 /dYN *^1 r dY V, 

^ mr “T“ ^ ^ 


k^ d^ Y 

“T 

Y dT 


i = 3, 5, 6, 7, 8 


a 4 ki) 

1 % 

dT^ 

- T "ar 

d k^ 

k^ E. 




2 

d^Y 

dl' 


RT 


7 J + 2 ( 


i — 3 f • • • fS 

d Yp 


d^Xf. d Yr, dX, 

^ . n. r ^ 


i =4,9 


■) (W") + 


”~7[r ^w~- 


\ 




dY. 


"W\ 


2 

d^X 


G 


dT 


[— g-, . ,^ . - .., 4 ,w + a] 

.(es+e^+'^+Sv- X, 

II —2 — 4- £r 4- W 4- a~t2 


T 


,2 ^ r ‘"s+^o+a-^v’ 

a 7 ^ H 


+ g + w + a] 
+ g + a] 


r %''’®c^‘^‘’’®v^ +g+w+a] 

I ■? 


4^0 2 


dT 




^(eg+o^+d+e^) + g+w+a^ 

j-_ _._ _ j 


f 9 w\ 2 


3 ^ w 2 

— _ _ — _ 

d T ^ ( e „+e^ +d+e ) +g+w+a 

[ o ] 



OHAPTER 4 


jRESULTS illl) EISCUSSIOIT 


With the objective function defined in Eq.(11), the 
optimal temperature profiles in batch reacoor are obtained 
for various reactor pressures using the combined method in 
which the control vector Iteration method is used for the 
first few iterations and then switching over to the second 
variation method till the change in the objective function from 
iteration to iteration is small. Eor starting the computations, 
initial guess temperature should be assumed. In the earlier 

studies it was found that ever, though the temperature history 

% 

is dependent on the initial guess tenperature assumed, the 
value of the objective function does not. Hence, we used an 
initial guess temperature of 500 ®K for the control vector 
iteration method and carried out all the computations. Eor the 
second variation method the near optimal temperature profile 
obtained from the control vector iteration method is used as 
the initial guess temperature profile. 

In the integration of the I)q.(20) or the Eiccati equation 
■ in the second variation method, with some values of the 
weighting parameters, we found that in some iteration, 
f coat t* where 0< t* < t^. At this point, as suggested 
by Bryson and Hot^ we stopped the integration of this equation 
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the objective function iomlnates and this would give smaller 
deviation of from This can happen only if the 

conversion of e increases. In lig, 3 , almost all the curves 

o 

indicate the use of high temperatures initially which must be 
brought do^m towards the end. Por low values of ^2* 
curves indicate the use of relatively low temperatures since 
the term corresponding to in the objective function has a 
very low value when compared to other terms and the effect of 
the other terms is to minimize the formation of side products 
and diethylene ^ycol content in the polymer. This can be 
accomplished by bringing down the reactor temperature. In 
particular, for a value of “2 = 0*5? the reactor temperature 
is seen to approach the lower limit set for the temperature. 

In I'ig. 4 j the increase in the of the polymer is plotted and 
it is seen that in all the curves, rises very fast and 
afterwards remains almost constant at an asymptotic value, which 
is attained in a relatively smaller time. It is implied from 
this behaviour of that in the optimization, reaches 
first and later on the minimization of the other terms in the 
objective function is taken care of. It is also found from 
lable 9 tint the value of is low for the low values of a 2 » 

In Pig, 5, for a value of «£ = optimal temperature 

profiles obtained by the combined method and the control vector 
iteration method are conpared. It is seen that the profiles 
are almost the same except in the initial and final stages 



EPPECT OP VARIATION OP a „ ON THE VARIOUS RESULTS 




on the optimal te 



lues of 






tf) ^-1000, ' too, pT=200fWiHg, Pnd *15 



Fiq. 6 CompQfison of th® jJn obtoined by the combined method and the control 
vector iterotion met 
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wherein the temperatures are slightly more or less, Ihis 
confirms our contention that the second variation method makes 
only small changes in the near optimal temperature profile 

obtained from the control vector iteration method. In Pig. 6 . 

•4 J.JC 

obtained for a value of "fehe ^ru-e- methods 

are compared and it is seen that the curves almost match in the 
initial stages, but the a^rmptotic values are slightly different. 
This is because of the chaises in the optimal temperature 
profiles in Pig, 5, 

The weighting parameter represents the relative 
importance of the diethylene glycol content in the polymer. By 
varying the optimal temperature profiles and the other 
results obtained are shown in Table 10 and Pig,?. In Table 10, 
it is seen that as increases, the conversion of e^ as well 
as the value of drops because of the domination of the first 
term in the objective function when compared to other terms. 

The result is that the moles of diethylene glycol in the polymer 
is minimized by having low conversion of groups. Since the 
conversion is low, the reactor temperature is also low as seen 
in Pig,?-, In this figure, it is also seen that the curves 
Indicate the use- of high temperatures initially for low values 
of but as increases, the initial temperatures begin to 

fall. In the final stages all the curves indicate the use of 
low temperatures. In particular, for a value of = 5000, 
the reactor temperature is seen to approach the lower limit 
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8 Cwnparison of the optimal temperature profiles obtained by the combined 
method and the control vector iteration method. 
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because the optimization problem reduces completely to the 
minimization of the first term in the objective function. 

In i'ig.8, the optimal temperature profiles obtained using the 
two methods are compared for a particular value of and the 
same behaviour is encountered as in Fig. 5. 

The vrei^ting parameter represents the relative impor- 
tance of the formation of side products in the reaction mass. 

By varying the results obtained are given in Table 11 and 
Big. 9. In Table 11, it is seen that as^- increases, the 
conversion of Eg and the value of both increases, but the 
increase in these values is less \jb.en compared to the corres- 
ponding values in Table 9, when ^2 varied. This shows that 
as is increased, the optimization problem reduces to the 
minimization of side products, thereby bringing down the 
conversion and to low values. 

5 

In Big ,9 as increases from 1 to 10 , the first three 
curves show that the initial temperatures ' are low, they 
increase to high ■values after sometime and then are brought 
down to low values towards the end. This behaviour is peculiar- 
to that expected since an increase in should lead to a 
reduction in the reactor temperature as more weight^e is given 
to the minimization of side products. This peculiar behaviour 
represents a complex interplay of the three terms in the 
objective function. 
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The effect of the pressure on the reactor performance 
has been studied and the results obtained are given in Table 12 
and Pig. 10, It is found from Table 12 that as the pressure is 
reduced, the conversion of E and the vd-ue of n increases. 
This behaviour is expected since the #M#-usion of the conden~ 
sation product is easier when the pressure is reduced, thereby 
favouring the polymerization reaction. It is found from Pig. 10 
that the initial temperatures are low, when the pressure is 
rsd.uced from 200 to 20 mm Hg, but for low pressure, the tempera- 
tures increase after sometime and then begin to fall towards the 
end. Shis behaviour is entirely different from that predicted 
by Eq.(4) in which the temperature is an increasing function 
of time. In Pig, 11 the optimal temperature profiles for the 
two methods are compared for a reactor pressure of 50 mmHg 
and it is seen that there is a considerable change in the two 
profiles. 

In Tables 13 and 14, the objective functions obtained from 
the control vector iteration' method and the second variation 
method are compared and ai can be seen that there is only small 
reduction in the value of the objective function obtained from 
the latter method. This is expected since the second variation 
method, derived by considering small perturbations in the 
value of the- state and control variables around the optimal 
profile obtained from the control vector iteration method ■ 
incorporates only small changes in the objective function 
obtained from the latter method. 
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11 Comjxirisonf of the optimal temperature profiles obtained by the combined 
method and the control vectoi iteration method. 
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li; can "be seen from Eq.(36) that the correction to he 
employed in the temperature profile depends upon the value 
of £ . It is found that the convergence is fast i.e., a 
maximum decrease in the objective function is obtained for a 
value of £ =1 for low pressures, irheroas for high pressures, 
convergence is attained for relatively low values of £ « 

All the results are generated using nBC-1090 computer and 
it takes about 30 minutes of computation time for obtaining the 
optimal temperature profile using the combined method. 



CONCLUSIONS 


TLe polycondensation stage of PEI formtion Las been 
assumed to include various side reactions in addition to the 
usual polymerization reaction. A simplified model for the 
flashing process is proposed in v;hich it is assumed that the 
flashing occurs at the end of small discrete time intervals. 
Luring this time, polymerization is assumed to occur in batch 
reactor. An objective function proposed earlier has been 
used in which there are three wei^ting parameters, 
and. . A suitable Hamiltonian is vjritten for the time 
Interval when the polymerization proceeds without flashing and 
the corresponding ad;io±nt equations have been written assumiing 
temperature and pressure as the control variables. V/hile 
computing the optimal temperature profiles, it v/as discovered 
no matter what the values of the parameters and 

are chosen, pressure always falls to the lowest limit, This 
leads to the simplification of the problem and in the following 
we have treated, temperature as the only ■ control variable 
for lower limits cf pressure. 

The optimal temperature profiles for a given lower limit 
of pressure were obtained using a combined method of the 
first and the second order techniques of opt5.mization. In the 
first order technique, the control vector iteration method is 
used till the change in the objective function is negligibly 
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small. It is well known that the convergence in first variation 
technique slows dom considerably as the optimality is 
approached. Because of this reason, the near optimal profiles 
obtained from the first variation technique is used as the 
initial guess for the subsequent second varia.tion technique. 

The solution of the Siccati equation in the second variation 
method poses both a memory storage as well as coirputational 

problem. The former has been overcome by storing the values 

_2 

at every time interval of 10 hr. For obtaining stable numerical 
solutions for the Hiccati equation, the integration has been 
carried out using a time interval of 2x10"*^ hr and obtaining 
the intermediate values of temperature, state, and adjoint 
variables by linear interpolation. In this integration 
sometimes, the Riccati variables, P, approached Infinity after 
some time of integration, Bryson and Ho suggests that is 
because of the existence of conjugate paths. Therefore, we 
stopped the integration at this point and assum©3 P to be zero 
for lesser times. This would lead to the correction of 
temperature profile for that portion of time, ^'ihere the values 
of P are finite. In this way, we could obtain stable 
numerical solutions . The weighting parameters involved in the 
objective function proposed for minimization, were found to 
have considerable influence on the optimal temperature profile 
we found that the combined method works better for obtaining 
optimal solutions rather than each method used separately. 



Computations suggest one must use Mgh temperatures initially 
which should be lov/ered subsequently. It has been suggested, 
that the high temperatures initially could be attained in 
practice by using PIT of high molecular weight as an inert. 

We also found that the optimal temperature profiles vrere found 
to depend significantly upon the reactor pressure and it was 
found that for low pressures, one must increase the temperature 
during the earlier portion of polymerization time which once 
again be lowered for large times to reduce the formation of side 
products. 
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1. COMPUTER PROGRAM FOR TBE GOMIROI VECTOR 
ITERATION METHOE 



4-* ****rt.**^* ***************** t***********^******************** 
□PTIHIZAITON OF PKT REACTOHS ♦ 

^ t- ¥ * * *****$ '¥*^*^¥**^* **^* *^* ********** ^^t**^'^****** 


PhPGPAP PSFS the CniiTROL VECTOR TTEkATIUN METHOD 


THIS 

vjPAniKMT *o:':r'nn[) FOK Fir^OING 
I'l A i'.\TCi< PLT 5?FArTuR 


THE OPTIMAL TEMPERATURE 


OH THE 

profiles 


Tits I'i/roni* IS ('ARRlEP OUT TILL THE CHANGE IH THE OBJECTIVE 
F^urTIOM r'of.n lTI,R-'vTlf*'i TO TTFRATTON IS SMALL 

IS'tTMf iruil. T£/4pFFATnHF IS USEO FUR THE IHTTTAL GUESS 

t ****** ****1f*****^* *'***** *****'¥*'¥*^**$3l(i*^********'*^^*^*** 


MAIN PROGRAM 

*it: t ****************************** ********************************* 


X STA'^h VARIABLES 

V Anoniw-f VARIABLES 

EO ACTIVATION ENERGY 

7:fi FREOUFMCY FACTOR 

Z RATE CONSTANTS FOP THE REACTIONS INVOLVED IN THE POLY- 

CnwOENSATlON STAGE OF PET FORMATION 
A «TI)TH of t'HE TIME IMfFRVAL USED IN INTEGRATIONS wITH A 

total TIME OF 2 HRS. 

EGO INITIAL concentration OF EG GROUP IN THE FEED 

PT TOTAL PRESSURE, atm. 

T REACTOR TEMPERATURE , DEGREE KELVIN 

T'lfcW TEMPERATURE AT 101 POINTS OF TIME 

TN temperature AT 2001 POINTS OF TIME 

K UNIVERSAL GAS CONSTANT 

ZP DERIVATIVE OF Z vv.R.T, T 

VDMT -molar volume OF DMT 

VG MOLAR VOLUME OF GLYCOL 

MOLAR VOLUME OF WATER 

VT TOTAL VOLUME OF THE REACTION MASS, LITRES 

Z1 Z DIVIDED BY VT 

HP DERIVATIVE OF HAMILTONIAN N.R.T, TEMP 

PI OBJECTIVE function 

VP DERIVATIVE OF VT W.R.T. T 

E epsilon-factor determining the MAGNITUDF of STEF TO BE 

TAKEN IN THE CORRECTION OF TEMPERATURE FROM THE PREVI- 
OUS ITERATUN 
EOPT OPTIMUM EPSILON 

AK VAPOUR PRESSURE, ATM 

ALPHAl - weighting PARAMETER FUR THE DESIRED CUNCE»TR- 
ATION OF DIETHYLENE GLYCOL IN THE POLYMER 
ALPHA2 - weighting PARAMETER FOR THE DESIRED VALUE Ot 
DPBAR OF THE POLYMER 

ALPHAl - WEIGHTING PARAMETER FOR THE FORMATION OF 
SIDE PRODTICTS IN THE PEACTUN MASS 
UPBAR NUMBER-AVERAGE MOLECULAR WEIGHT OF THE POLYMER 

**********************,******************************************** 


DIMENSION X(7, 2001), V(4. 2001), Z0(9),E0C9) 

DIMENSION TC20).TNE«(10i),TRCl01),TNC2001) ^ ^ ^ 

DIMENSION ZC9,2001) ,ZPC9,101),HP(101),EP(1S),QI(15) 
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DXyE.iStiiii ,?n»tT(2ooi),vc;(200lJ,VW(2u01),7T<200U 

l-:3(1U),FiClO) .21(9,7001), nPEAR(2001),OF(l5j 


Cf'..!Mfjf;/A ! /3CA(h'i/h2/A/^^/Z 

C -j fc / A 4 / V v; , V i ) ' * T , y t; , V l' / h 1 / EU , 20 , R 

cn,., iiF!/3:,//,i /‘-p/TM/n?/^ 

cn-' <'j:.vah/zp,hp 

C'’-(MU*i/AO/FGO/A lU/TfJt:W 

Cl’ /oT,H 

Cfi ■riiVi/Ar/nloTF 

Cnv,M0'!/‘'il2/inO'AF 

PC/J, -M* 


K>';cc=i' 

.-;(;ii=i,r, 

,= -}.!> 

!YCF. 540 

FOIMaTI 1 X, 'GIVE lOE VALUES OK ALPHAl , ALPHA2 , ALPHAS IM D.LINEti'J 
KEAH (5,*) ,A!FHA1 ,ALPtlA2, alphas 
M nTE(t}0,3) ,ALPHA1 .ALPhA2, ALPHAS 

FORM at ( I X, 'ALPHA1 = * ,F15.5,2X, 'ALPHA2-',F15,5,2X, 'ALPHAS^' ,f Ib.b) 
0=1 Cl 
TYPF 67 S 

FfU'OATdX, 'GIVE TOTAL PPESSHHF') 

UFAn(S,*),PT 

PT=0.0 

A=1.0/10U0.0 

K=l,997 


140=0 


KEAr)t24r*) ,(TM(I),l = l4ll) 

WPlTECbS.I l99),TNCn,PT 

FOR MATO X, 'GUESS TEMP: = ' ,F1 0 . S , 'TOTAL PRESSURES FI 5 , U •» 

UPErl(UNIT=23) 

rHEWCU=TNCl) 

Jl=2 

DO 444 J=l,10 
D=1,0 

DO 444 1=1, to 

TCl) = (TNtJ + l )-TN(J))*D/10,()+TNCJ) 

TMEW(U1)=T(J) 

U=D+1,0 

Jl=Jl+l 

CONTINUE 

£0(3)=lftb00.00 

EO(4)=290Ou,on 

E0(5)=29a00,00 

E0(6)=29800.00 

En(7) = 17t>00.00 

£0(85=17600,00 

EOC9)=37BOU.OO 

ZO(3)=13.4i 

20(4)=17.54 

Z0(5)=17,54 

Z0(6)=17,S4 

Z0(7)=i3.855 

Z0(8)=13.R55 

Z0(9)=22,004 

U0=M0+1 

wRlTE(60,16) 

WRITEC60,12),N0 
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FC!P-MAT(IX,120<’1H-) J 

TYPE 12, an 

{•i1HYaT( 1 A., 'NU=' ,12) 

■/iRlTECf.n, lb) 

-.PJThCbi-, lb) 

(’AbCi'LAT LfPJ (IF MOfjAR VULUMES AMD RATE CONSTANTS AT 2001 

P^-Uf'i'S '■'F rTf.ip: 

iH> 1b i = i,S» 

Z(1 ,1 )s:6O.0*HXpr-Fu(I)/fH»TMtWtU)+7aCl)) 

;.PCT, 1 )=r;rKl I. , 1 VfF^'iEbWtl )*TNEW(1 ) ) 

C”! TINliF, 
r;. (I ) = 1 1 1 

T'{ 1 ) = iMi,b*C1 .0 + U.0i)l4*(TN( 1)-413.0>)/1000.0 
0 1 1 j =t. 0 , f- »• ( t , (1 + 0 . on 1 A f Tf j ( 1 ) -4 1 3 . 0 ) ) / 1 000 . 0 

Vi. ( 1 ) = ( 1 p . 4 22 + 0 .025 * (Tn f I )-41 3 .0) ) /I 000.0 

j ^ — 'I 

3=1,100 

n=i .0 

Dit 14 3 = 1,20 

r C 1 3 = ( Ti'ii'.F ( J+ 1 )-TNEW c J ) 3 *D/20 .0 + INEWC J3 
u(l 15 I1=3,P 

X(3 1 ,.]l3=no.U*EXPC"EUCll)/CR*TCI)) + ZO(ll)3 

COi.TlNijF 


roUDsTCTJ 

vnfiTCJl 3 = 101, 5*(1,0 + 0.0014*(TNCJ13-413. 0)3/1000.0 
/G(313=OO.b*C 1.0+5.0014* (TN (31 3-41 3,0) 3/1000,0 
3 = ft 5. 422 + 0. 025* (IN (31 3 -41 3. 03 5 /1 000.0 
3=r» + 1.0 
31=31+1 
CnNTlNlJE 

IMlTIAi, CONDITIONS CJN TFE STATE VARIABLES 

X(l, 13 = 1.0 

X( 2,1 3=0.0 

X( 3, 13=0.0 

X(4,l)=0.0 

X( 5, 13=0.0 

X(b, 13=0,0 

XC7, 13=0.0 

CALL FUTJcaiA) 

NRlTC(bO,67P3 , ((X(I, 33,3=1 ,2001,1003,1=1,73 
F0RHAT(2 a,'X is; ’,21F10.3) 

WR1TE(0{^,16) 

TYFF 791,N(j,DI5TF 

FORMATdX, 'FUR N0= ', 12 AMOUNT FLASHED= ' , Fl5 . R 3 

wRITE(b3,7Ql) ,Nn,DISTF 

1>ISTF=0.0 

FINAL C0r4DITinNS ON THE ADJOINT VARIABLES 

V(l, 20013 = 0.0 

V(2, 2001 3=0.0 

V(3,2001)=ALPHAl 

V(4, 20013=0,0 

MP=(X(1 ,20013+X(2,2001 3+X(3,2001 3+X(4,2001)3/2.0 
VT(2001)=VDMTC2001 3*MP+VG(2OO13*X(5,200l3+VW(200l)*X(6,2O01) 
Z1(4,2001)=Z(4,2C01) 

Zl(9,2o0l)=Z(9,2001) 

DO 191 1=3,9 

IF((i eQ. 43,UP,(I.EO,933GOTO 191 
Zt (I|2001)=?Z(l,200l5/VT(200i3) 

CALL RUNG2(ALPHA2, ALPHAS, A) 
bRlT£(&3,6003,((V(l,J),J=t ,2001,1003,1=1,43 
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wPITEft.i, to) 

TYPH 12 ,l '0 

.) 0 (! FnHfiAT(2X, *V Is: •,?1E11,3) 

PPI=PT 
CAl.l, SiNFSit 

PT=ALPH,Altxf 3,2001 ) + ( ALPtiA 2 / 225 , 0 ) ♦SSX ( I ) + ALPHA 3 » (SSXC 2 ) 

1+S.SXC4)) 

iJFCf3U)=Pl 
TYPE 23 , FT 
rtRlTKC60,?B),Pl 

2 B FOhMATClX, 'PI=',F 20.51 

type U3 
RFAP( 5 , 4 b!S) ,C 1 
IF^Cl.FO. ';•)* jFilTlJ 42 
IFt'i.l.E J. nFfJTiJ 1510 

TF((ABS(UF(fJtl-l)-nF(Nn))),LE,O.OOUGUTU 42 
1516 IFCHti.GT.g) GHTn 42 

FO 32 ,100 

DO 30 1 = 1,7 
K( 1 ,J + 1 ) = X(I,J* 2 ()+i ) 

iFtr.GT. 4 )«nTn 30 
v(i,j+n=v(i,j» 2 o+i) 

30 COUTIHUF 

VT(J + 1 ) = VT 1 J 1 ‘ 20 + 1 ) 

32 CONTINUE 

DO 33 1 = 3,9 

IF((I.E 0 . 43 .r)P,(I,EO, 9 ))GOTO 33 
1 P=(X( l.n+Xf 2 . 1 )+XC 3 ,l)+X( 4 ,l))/ 2.0 

VPs( ( (MP 3 * J 91 , 5 » 0 . 001 4 +X( 5 ,l) » 60 , 6 ^ 0 , 001 4 +X( 6 , 1 ) * 0,025 3 
I/IOOG.O) 

ZP(T, 1 3 =(ZPf i,l 3 /vr(l) 3 -((Z(I,l)/(YT(l 3 *VT(l)y 3 *VP) 

33 continue 

DO 31 0 = 1,100 
UO 31 1=3 9 

Z(I,J+ 1 )=^U, 0 *EXPC-EO(I)/CR*TNEW(J+ 13 )+ZUCI 3 ) ^ 

ZP(I,J + l 3 =En(l 3 *ZCIrJ + U/CH*TNEWlJ+ 13 *TNEWCJ+ 13 ) 

IF((i EO. 43 ,OP.(l.E 0 . 93 )GOTn 31 , ^ 

MP=(X(l,J + 13 +X( 2 ,J + 13 +X( 3 ,J + l)+X( 4 ,J+U 3 / 2,0 

VP=( ((MP 3 * 191 , 5 * 0 . 0014 +X( 5 ,J+ 1 )* 60 , 6 * 0 . 0014 +X( 6 ,J+ 1 )* 0,025 3 

1/1000.0 j , 

ZPCI ,J + 13 = (ZP(I, J+1)/VT(J + 1 3 3-C(Z(I,J+l)/(VT(J+l)*VT(J+l 3) )*VP) 

31 COfJTINUF 
CAGE DHU(M 3 

TYPE 22 . (HP(J 3 , J =1 ,M) . ^ 

22 FORMATClx, 'THE VALUE OF HP IS ' , t 01 El 5 , 4 3 

WR 1 TK( 60 . 22 ),( HP ( J 3 ,J= 1 ,M) 

T Y P E! 1 2 fi 0 

18 FORMATcl X, 'THF STARTI NG EPSILOH= ' , FI 5 . 4 , ' IN THE INTERVAl-s'.FlS -4 

121 TYPE 29 

UO 601 K 2 =l ,10 
EPCK 2 )= 0.0 
QI(K 2 )= 0.0 
501 CONTINUE 

Kl = l 

29 FORMATClx, 'OIVE THE VALUE OF STARTING EPSILON AND INCREMENT . *3 

READ( 5 ,* 3 ,E,FNEW 
E 3 CN 0 )=E 
F 3 CNO>=PHE '4 
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.'i«lTI-:Cf>f;,700i) ,F3(Nn) ,F3CN0) 

7003 FHRMiVTC f X, 'THE ST^lt^TlMG FFSILO!J=' ,Ft5,5 , 'THE INCREHEHT=' .F15-5 » 
lOf ? Fi)HMAT(FiO.IO) 

57 un 51 J=1,101 

TR {..I ) =T?^KW i H ) -E*HP ( J ) 

IFCTHCJ) .T,T.bO0,0)GOTn 5 
TRCJJ=bOo,o 

5 lEfTHC J) .€T, 400.0) GOTO 51 

TR(. Jj=400,o 
51 CHf^TlNtir 

fin 00 1=3,9 

^ T,l)= 00 ,o*EXP(-EnCl)/(H*TRCl))TZOCI)) 

bO CUiiTlNUE 

7r>‘'T(U = iQi ,5* ( 1 .0 + 0.00t4*(TH(l )-413.0) )/100U,O 
tfG(l )=o0.6»f 1.0 + 0.0014 ♦(TRtl )-41 3.0)) /1 000.0 
i/W( t ) = CIO, 422 + y,0?51^CTK(l >-413,0 ))/1000,0 
J 1 = 2 

TMUJ=TRCn 


53 


52 


166 


B2 


72 

123 


on 52 J=l,100 
11=1,0 . 

on ^2 1=1,20 

T(I)=fTP(J+l )-TR(J3)*n/20,0+TP(J) 
on 53 11=3,9 

ZCn jJl)=6U.O»EXPC-EO(Il)/tR1‘TCI))+ZO(Il )> 
ffJCJi)=TCT) 

VDHT(vll) = 191,5*( 1,0+0,001 41' CTNCJ1)-41 3. 0))/tyo0.0 
VC( Jl)=60,b*(l .0+0.001 4* CTNCJl)-4l 3,0)) /lOOO.O 
VW (J1) = C 19. 422 + 0.025* (TN (.71 )-4l 3.0 ))/1000.0 

n=D+1.0 

Jl=Jl+t 

cohtinue 

CAbL RUNGl(A) 

DTSTF=0.0 
TYPE lb6,Kl 

FORMAT (IX, 'NO. INSIDE THE LU0P=',12) 

PI I=AU'riAf*y( 3,2001 ) + CALPHA2/225,0)1'SSXCI)+ALPHA31‘CSSXC2) 

1+SSXC4)) 

EP(Ki)=E 
TYpF 82 

format ( tX, 'MORE?') 

READ(5,455) ,B1 
IF(B1,FO.'H*)GOTO 72 


QT(KI)=P1I 
TYPE 7 


C*********** 


YPE 71,EP{K1) ,0ICK1) 

TYPE 82 

READ(5,455) ,Bt 
IFCPt.EO. 'f.')GuTO 72 
E=E+FNEW 
El=Kl+l 

IF(K1,LT,1 UGOTO 57 
TYPE 12, NO 
TYPE 123 

FORMATCIX, 'DO YOU WANT TO GO BACK IN E .OPT.^OOPf Y/N > ?! ' J 
HEApC5,455),Hl ^ ^ 

♦:^***l******!******c*;|;*^:J***)(:f*!!***********!** 
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b 5 


wFnTfc;CbO, 1 b) 

TYPK 

•(RITtr’J.bO^bS) 

TYPF 65 

KnHMAt(15X, 'I' ,10X, 't;' ) 
pn 7<) = 

/jfUTfc;fbO,7n,otfKn ,EP(Ki ) 
<91TEC5,71 },ytCKl),HP(Kl) 
F(.iH 6 ATt 7 X,Pr;. 4 ,i; 


71 F(.tH6ATt7X,Pi7.4,iX;‘=‘i7.7) 

7 0 C(H''Tifnjp': 

61 FnR‘1AT(2X, 'FUPT=:' ,E20.10J 

TYt’f- 1 00 b 

lOOB PnOflATdX, "GIVE THE YAl.t»E HF FHPT') 

KFAn{5,=l'),rtL 
1u.j9 FORMATcFI 5,4) 

CnpTsAL 

4RITLCbO,lD) 

WblTECbO-64) .FOri 

C THIS PART CORPErTS THE OOP TEMP. PROFILE USING EOPT 

on 40 d=l,M 

TflhW(J)=TNEWCJ)-EOPT*HPCJ) 

IFCTHFwC J) .LT. 600.0)GOTn 17 
TMEW(J)=600,0 

17 IFfTHEWCJ) .GT.400,0)GOTn 40 

TNEWCJ)=4U0,0 

40 cnoTiNut: 

c f^*’^*****t* ********** ** 1 ^* 11 *^* **^***!lli^>f*****^ *********** 

41 FORMAT! IX, 12) 

44 FfJHHATClX,101F11.6) 

TYPE 45, lTr4FW(J) ,J=:1.M) 

45 FORMATCJX, 'TNFH IS; *,101F?,2) 

APlTECbO,45) , (TL'EWCJ) ,J=1 ,M) 

113 FOKMATdX, 'DO YOU ftArtT TO PROCEEDC Y/N) 1 ' ) 

455 FORMATCAl) 

TYPE 113 
READ(5,455).F1 
1F(F1,E0. 'Y*)GDT0 499 
GOTO 42 

499 CONTINUE 

49 GOTO 11 

42 WRITEC23,549) , (TNEWCJ) ,J=l ,M) 

00 46 .J=1,M 

IFCTNEWCJ) .LT.bOO.fDGOTO 4 
TNE:W(J)=600.0 

4 1F(TNEW( J) ,GT.400.0)GOTO 46 

TNEW(J)=400.0 

46 CONTINUE , ^ 

TYPE 549, (TNEW(J). 0=1,101) 

549 FOBMAT(1X,101F7.2) 

646 F0HMATaX,F7,2) 

DPBAR(206i )=(EG a/CX(l ,2001 )+X(2,2001 )+X(3,2001)+XC4,200l O ) 
WRITE ( 59 , 899 ),CX(l,J),J=lr 2001, 20) 

WRITEC59,16) , . 

WRITE (59, 899), (X (2, J), 0=1,2001 ,20) 
wRlTEC59,J6) 

WRITEC59,899),CX(3,J),J=1,2001,20) 

WRITEC59,899) ,(X(4,J) ,0=1 , 2001 ,20) 

WRITEC59,16) 
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•■IPITLCS^', Ih) ' ' ' 

/FRI r^-:(59, Ryq), (X (7, J),J=i, 2001 ,20) 
rti'!ITE(b9,lb) 

■/riTKCbO.syq), (r<PRAR(I),I = l,2001 ,20) 

VJRXTE(5R, lb) 

8«9 FUHH4T(1X,'}« is: ',101E15,6) 

w!>ITECb0,12),N0 
STdP 
END 

****** **^*************t****^:********4***********V*** 
THIS SHPHnUTINE INTEGRATES ^HF STATE VARIABLE EOUATTQNS 
USING FnuPTH-PHPER PUNGF-KUTTA METHnD 
SUriROUTiNF FUNGI (.A) 
niMFMSTPN X(7,200n,XORGC7),XT(7) 

DIMENSION XnRGl(7),XR(7) ,TN(2001 ) ,En(9),ZO(R) 

DIMENSION VDMT(200I),VG(200t},VW(2001),VTf200tJ 
DIMENSION ZI (9,2001) ,ZIC9),DPBAR(2001),ZC9-2001) 

DIMENSION CP1(7) ,CP2(7J ,CPS(7),CP6(7),CPC7LfPC7) 


C 


99 


12 

5555 


101 


1111 


COMMnN/A2/X/A3/Z/A4/VW,VOMT,VG,VT/ei/EO,ZO,R 

C0MM0N/A5/Z1/A6/TN/RUNGU/FP,XK 

COMMUN/Rl /XI,X0PG 

COMMUN/A9/EGO 

CnMMuM/FFl/PT,H 

COMMON/Bb/DTSTF 

CPMMnW/AlZ/DPBAP 


REAL MP 
REAL Mpp 
A=A/5.n 

DO 100 Jl=l,2000 
J=J1 

MP=(X(1 ,J)+XI2,J)+XC3,J)+X(4, J))/2.0 
DPbAR(d)=( (EGn/NP)/2.0) 

VTCJJsVDMTC J)»MP+VG(J)*X(5, J) + VW(.JJ*X(6, J) 
DO 99 9=1,7 
XORGl (N)=X(?^,J) 

CONTINUE 

N=1 

D=0.0 

DO 12 1=3,9 
ZUI)=Z(1,J) 

continue 

CONTINUE 

DO 101 N=l,7 

XnRG(N)=X(N,J) 

XR(N)=X(N,J) 

continue 

Tl=TN(J)+D*(lN(J+i)-TH(J)) 

Z(ll jUeoTo^EXPC-EuCD/CR+TU + ZnCI)) 

continue 

xi=191 ,5*CI,0+0.0014*CT1-413,0))/1000.0 
X2=60,6*(1.5+0.0014*(T1-413.0))/1000,0 
X3=(1§,422+0.025»(T1-413.0))/1000,U 
MPP=£X6 rGC l)+X0RG(2)+X0RG(3)+X0RG?4))/2,0 
VTO=Xl*MPP+X2^XC5,J)tX3*X(6,J) 

DO 2222 1=3,9 



m 


2222 

1U2 

1 Q 3 

1 0 6 

10^ 

999 

578 

34982 

107 

100 

C 

c ■ 


7222 

^1 ( T,»n=zr i,in/vTn 

COfiXlNUF 
4l C4,vJ)=Z{4,U) 

Z1 (9,.J) = ZC9,U) 
l>=I)+0.2 
CAl,L FUTNl tJ) 
on 1(.>2 1 = 1,7 
CFK n = A4.FPfI> 

JC9tI) = A(JHn(I)+0.5tCPiCl3 

CnflTINIlF. 

CALL EOTNUtT) 

UO 103 1=1,7 
CP2(X3=A»FP(I3 

CALL FQT.’Jl ( J) 
on 106 1=1,7 
CP3C1)=A»FP{I) 

XR(T)=xni«:::(T)+CP3CI) 

COfTIOUF 

CALL EUT61(J) 

DO 105 1=1,7 
CP6( t)=A*FP(n 

COLTINOF 

■1 = 0 4-1 

IF(N.L6,5) COTn 5555 
on 999 1=1 ,7 
X(I,J+1)=X(I,J) 

X(I,J1=XURG1 (I) 

CONTINUF 
DO 578 1=3,9 
Z(I,j3=ZIfI3 
COOTiNUF 
T=TNCJ) 

IFCXCS.J+n.GT.O.O) GOrO 34982 

Al = Zl(3,J)*XOPGCl)4XORGCl) + (Zl (7 , J) ♦XORGC I ) ♦XORG C6 ) /2 . 5J 
B=Z1 (3,J)=*'4.04CFGO-XORGC1)-XOPGC2)-XQRG(4)) + 2.0*ZU5,J)»XOHG11 
l+2.0*Zl (7,^J»xnRG£l) 

C1=CA1/B) 

C2=FXP((-B)*An 

XC5,J+l)=ri*(l ,0-C23+X0RGC5)4C2 

CONTINUE 

DO 107 1=1,7 

XORr,(T)=X(I,iJ + l> 

CONTINUF 

call flashcdist,nox,t,j) 

IF(NOX.Ea.O) GOTO lOO 
X(5, J + l )=X1(5) 

X(6,.H-n=XI(b) 

X(7,J+n = XI(7) 

DISTF=DISTF+DIST 

continue 

A=A»5.0 

RETURN 

KNf) 

THIS SUBROUTINE CONTAINS THE STATE VARIABLE EQUATIONS 
SUBROUTINE EQTNl(J) 
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c. 

c. 


CUmMUN/AS/ZI/RUHGI 1/FP,XR 
CnMM.o?^/A 9 /Ef:a 


c 

c. 

c 

c 


55 


12 


i'P(n=-2.0*7H3,J)*(XR(l)=l'XR(l)-(C(2.0=t<XRf51>70.5)*(l,0-XRll J-XR 
n2)-XP(4))p-z|(4,J)*XF(U-2.0»Zl(5,J)*XRCl J*XR(5)-2.0»7-J 
lXPa)*XP(l)+ZH7,J)»{2,0*XRC2)»XRC5)-(XR(l)+XR(t>)/2.5Jl, ^ ^ 

1-Zl C8,J)*tXRC2)4=XR(l)-(XR(6)*(RGO-XR(n-XR(2)-XK(4))/1.25i) 
l-Zt (3,J)*XP(4)=t=XR(l) 

FP(p=Zl (4,J)fXR(n+ZiC5,J)»2.0»XR(l)*XR(5) + Zl(6.J)*XRCn»XR 

l(n-Zl( 7 .J)*(t 2 . 0 *XR( 2 )»XR( 5 ))-CXRa)*XR( 6 )/ 2 , 55 ) 

l-Zl C8,Ji%(XRf2)*XRCl3-((XR(b)/t .251*fRGO-XR(l)-XRt2)-XHJ4l)i) 
1+ZU9, J)*(FGO-XRf 1 )-XR(2)-XRf 4))/2.0 
FP(3)=2,y*ZU5.J)fXR(l )*XR(5) + Zl(6,J)»XRU)*XR(l ) 
FP(4)=Zl(9,J)*tFGO-XR<l')-XR(2)-XRr4))/2,0-Zl(3,JJ*XRC4J»APCl i 
FP15)=:Z1 (3, J)»lXR(i)*XR(l )-4,0*XK(5)^CEGU-XR( t J-XRt2)-XRf4)) J-Zl 
l(7,0)»C2.0*XR(2)^XR(5)-CXR(l)*XR(6)/2,5))-2.O»Zl<5#J3»Xk(i)*XHf. 
15) 

FPC6)=Z1 t7,J)*(2.0»XR(2)*XR(5)-(XR(l)*XR(6)/2,5>)+ZlC8,J)*tXR(2) 
l*XRCl)-tXR(b)/1 .25)^(EG0-XRC13-XR12)-XR(4))) 
FP(7)=Zl(4,v1)*XRCl}+ZU3,J)»XRCl)*XRC4) 

RFTHHN 

ERi) 


THIS SUBpnUTINE IS USEU FOR SDEVING THE ADJOINT VARIABLE 
EOUATIUNS USING FOURTH-ORDER PUNGR-KUTTA METHOD 
SflBRQUTlNF RUNG2(At.PHA2, ALPHAS. A) 

UTMFNSIOK XC7.200U,V(4.2001),ZC9.20U1),Z1(9,2001) ,ZIC9J 
dimension vnHG(4) ,VnRGl(4),FP(4),VRC4),EOC9),ZU(9) 
dimension CPC4).CP1(4).CP2(4),CP3C4).CP4C4) 

DIMENSION VWC206u.VOMT(2001).VGC2OO1),VTC2O01) 

DIMENSION VAPf73,TNC2001),TNEKC101),XR(7) . 

’COMMfJN7A27x7Al7z ' — — ^ ^ 

COMMON/ A4/i/W.VDMT,VG,VT 
COMMON/ A5/Z1/A6/TN/A7/V 
COMMON/A 30/TNEW/FFl/PT,H 
CaMMON/A9/EGO 

C0MM0N/RUNG21/FP.VR/E1/E0.Z0,P 

REAL Mp 
REAL MPP 
A=A/5.0 

DO 200 JI=l,2000 
J1=2002-JI 

MPsCXCl .J)+X(2,J)+X(3.J)+XC4. J))/2.0 
VTCJ)=VDMT(J)*MP+VG(j5»Xt5, J)+VW( J5*X(6,J) 

DO 55 1=1,4 
VORGl(I)=V(I,J) 

CONTINUE 

Nsl 

U=0.0 

DO 12 1=3,9 
ZI(I)=ZCI,J) 

CONTINUE 
DO 399 1=1,7 



99 


5S55 

I n 3 

3 

a88B 


COivTIfJUK 

^n^ tinuh 

tm 1111 i=i,9 

z f t » >,' ) =flO . 0* F.XP C -Efi ( I ) / ( R* T ) fZU( 1 ) ) 

CffjTIWljF- 

un 3yp 1=1,7 

^ c 1 4 j ) f;X C T , J ) -u* ( X ( T , J ) -X ( I , vJ- 1 ) ) 
COfjT 1 h UK 


a 1 = 191 .5* f 1,0 + 0,001 4»{i’-4t 3,0) j/iuuo.O 
X?=60.b»(1 ,0+0.001 4*{T-41 3.0))/1000,0 
X3=(19.422 + 0.025 + fT-413.0))/10C>0.0 
riPP=(XCI ,J)+XC2,J)+X(3,Jj + X(4f J))/2,0 
VTO=Xl+MPP+X2*X(b,J)+X3*X(b,J) 
m B«RH 1=3,9 

rK( {I.F.0.4) ,OP. CI,EQ,9))G0T0 8888 

Z(l,J)=ZtI,a)/VTO 

CnOTINliK 


l.>=D+0,2 
UO 201 N=l,4 
/nRn(N)=vcr!,j) 
v}u^^)=v(^,J5 

201 CCU.TINUF 

CAI.L KgTf^?(ALPHA2,ALPHA3,J,T) 

DO 202 M=l,4 
CPI (N)=A»FP(N) 

VP(rO=vnkGCN)-0,5*CPl{N) 

2<'2 continue 

CALL E’yTN?(ALPHA2,ALPHA3,J,T) 

DO 203 M=1,4 
CP2(N)=A*FP(N) 

VR(M)=VnRG(M)-0,5*CP2(N) 

203 continue 

CALL EoTN2CA1,PHA2,ALPHA3,J,T) 

DO 206 N=l,4 
CP3(M)=A+FP(N) 

VRCN)=vnRG(N)-CP3(N> 

206 CONTINUE 

CALL EQTM2CALPHA2, ALPHAS, J,T) 

DO 205 N=l,4 
CP4CM)=A+FP(N) 

CP(N)=(CP1 (N)+2.G*(CP2CN)+CP3CN) )+CP4CN) 1/6,0 
VCN, J)=V0PG(N)-CPCN) 

205 CONTINUF 

N=N + 1 

IF(N,LE,5)G0T0 5555 
DO 56 N=l,4 
/(W, J-n=VtN,J 
VCN, J)=VUPG1 CN 
56 continue 

DO 578 1=3,9 
Z(1,J)=Z1(1) 

57U CONTINUF 

1)0 397 1 = 1,7 
XCl,iJ)=XB(t) 

397 CONTINUE 

200 CONTINUE 

A=A*5.0 
RETURN 
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ADJOiHT EC'UATtQNS 
fOTh2tALPHA2, ALPHAS, J,TJ 
Xf?, 2001), Z(9, 2001), FP(4),^/R(4) 

Uli-E^iSinr! TL( 200 n,TNEw(' 10 l),VA'>C 7 ) 

CnMMiiN/A 2 /X/A 3 / 7 . » V ^ 

cc '^MON/ROnG 21 /FP,Vk 
C 0 i-,Mnfi/Ab/TM/A 30 /TNEW/FFl/PT,H 

f*F THE FIRST derivatives OF “F” W.R.T."X*' 

I ^ 1 1 -I C I J } 

VAP(5) = 1 (1.0** (21 .61-3720,0/T-4.042''=ALOGIO(T))/760,0 
<^8. 064 103-1 757,953/ (T-33, 274)) /760.0 
P=PT/(2.0*1VAP(6)-PT)) 

/JXl=P? WX2 = PrwX3=:P?rtX4 = P 
U=PT/(2.0*CVAP(5)-PT)) 

GXl = g?r,X2=:0;GX3=9;GX4=Q 

pFlXl = -4.0l7.(3, J)*X(i .J) + ( (4.C*Z(3,.J)/0,5)*(-XC5,J) + (EGU-Xtl . J) 
|-X(2,J)-XC4,J))*GX1 ))-Z(4,J)-Ct2,6*Z(5. J))=r(X(5,J)+XCl »u)+GXl 
U)-;Up*XCfa,J)*X(l ,J)+2,0*4(7,J)*Xt2, J)*GXi-C(Z(7,J)/2,p J*txtl . 

I J)*WX1+X(6,J)))-Z(8,J)=»'X(2,u 5 + ( tZ(8,J)/l,25)t(-X(6,J)tlFGiJ-A( 

II ,.T)-X(2,J)-X(4,J))*WXl))-Z(3,J)»X(4,a) 

pf;lX2=a4,0*ZC3,J)/0,5)*l-A(5,d) + (EGO-XCl,J)-XC2,J)-X(4,.i) lfGX2 
l))-2.0*Z?5,J)*^(l,vJ)*GX2 + CC2.0*Z(7,J))*CXC5,J)+X(2,J>’fGX2))- 
1 CZf7,J)+X(l ,J)/2.5)*WX2-ZC8-J)*XCl ,J) + C(ZC8lj)/1.25)<'(-X(b,y) 

UtEGO-XCI ,j5-X(2,d)-Xt4,J))iwX2)) 

1)F1X3=0,0 


l*WX4))-Z(3,a)*XCl,J) 

L)F2X1=Z(4, jHC (2.0*7,C5, J) ) ♦CXC 1, J) *r,Xl +XC5, J) ) ) + 2. 0*Z(6,d) ♦Xtl .J 
l)-2.0*Zt7,J)*X(2,J)*GXl+(CZC7,J)/2,5>*CX(l,J)*WXi+X(6,Ji) 

1 J)*X(2,J) + ( (Z(a^J)/l,25)*(-Xf6,J) + fEGU-X(l,J)-XC2,J)-Xl4,.i) J* 

.iF2X^=2!o*^.(5^j5*XCl ,J)*GX2-C(2.0 + Z(7,J))*(X(2,J)*GX2+Xlb,d) J) + 
1CZ(7,J)*X(1 ,J)/2,5)*WX2-Z(8,J)*X(1 ,J)+CCZt8,J)/l,25)*l-X{h,d)+ 

1 (EGO-XCl , J)-XC2,J)-XC4,J))*WX2))-Z(9,J)/2.0 
UF2X3=:0,0 

DF2X4=2.0*Z(5,J)*X(l ,J)*GX4-2,0*ZC7,d)*X(2,J)*GX4+(Z(7,Jl*X(l,J) 
l/2,S)*WX4+( (Z(8.J)/{.25)*t-X(6,J)+(EGU-X(l,J)-X(2,J)-Xl4,a))i 

lWX4))-ZC9.0)/2.0 

DF3Xl=((2,0*Z{5,vn)*CXt5,J)+XCl ,J)*GXl))+2.0*Z(6,J)*X(l».O 

DF3X2=2.0*ZC5,J)*X(1 ,J)*Gy2 
1)F3X3=0,0 

UF3X4=2.0*Z(5,J)*X(1 ,J)*GX4 
DF4X1=-Z(9,J)/2,0-Z(3,J)*X(4, J3 
0F4X2=-Zt9,J)/2.0 
[)F4X3=0,0 

OF4X4=-Z(9,J)/2.0-ZC3,J)*X(l,J) ^ 

“CALCULATTofToF-’W^ ' ' ^ ^ ^ ^ - 

FPC13=(2.04Al.PHA2*EGU*((EG0/X(l,J))-i5.0)3/(225.0*X(i,J)*XCl,a)) 

1-(DF1X1*VP( 1>+DF2X1 ♦VR{2)+DF3X1*VR(33+DF4X1*VRC4)) 

FP{2)=-2.0^ALPHA3*X(2,J)-(DFlX2*VR(n+DF2X2*VR(2)+DF3X2»VRt3) 

1+UF4X2*VP(4) ) , 

FP(3J=-(UFlX3*VR(n+DF2X3*VR(2)+DF3X3*VR(3)+DF4X3»VRt4)) 

FP(4)=-2.0*ALPHA3*X(4,J)-CDF1X4*VRC1)+DF2X4»VR(2)+DF3X4*«RC3) 

1+DF4X4*VPC4)) 

RETURN 

END % 



C 
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c 

c 


30 


c 


c 


c 


1008 

99 


n SUB<fHTI0h C/lliCULATES THF! DKRIVATXVE OF THE HAMILI- 

(Jf'lA)'! '^',k,r. TEMP 

5n(.>RQUTlr4E 

HPdOl ),X(7,2001),ZPC9,10U,V(4.2001 1,2(9,2001 1 
OTHii^Sinw VAF(7'»,ThiEw(iOl),VT(206l),\fW(206n,VGC20Ol),VOMTt2U013 
CO-.On^'/A8/ZP,fip/A2/X/A7/V 

Cnr’,MUf4/A9/EC:0/FFl/PT,H/A30/TNFW/A3/Z/A4/VW,VDMT, VG,VT 
UO 30 
0=01 

T=TMEVi(J) 

'(APC5) = 10.0=)‘=K21.61-3729.0/T-4.042»ALQG10(T))/760.U 
\/AP(b 8. 06410 3-1757, 853/(T-3 3,274))/760,0 

A=VAf’f5)l'C ( 3729.0/ tT»T))-4.042/T) 
ri=\/AP(b)*( 1757.853/(1-33,2 741=»‘»2 3 
C=-(X(5, J)*A+X(6,J)*B) 

GT=C/(VAP(5)-PT) 

ivT=C/CVAP(t>)-PT3 

iiP(J)=0.0 

HP(J) = V(1 ,J) + C-2.0*ZP(3(J)=<'(X(1,J>=1=X(1 ,J)-((4,0^X(5,J)/O.5O)J(l, 
lU-XC 1 ,iI)-X(2,J3-X(4,J))l)-ZP(4,J)»X(! ,J)-2,0*ZP(5,J)»XCl ,u>fX(5 
1 ,J)-2,0*ZP(6.J)»X(1 ,J>^X(1 ,J3+ZP(7,J)4(2.0*X(2,J)*XC5,O J-lXti .J 
l)tX(6,J)/2,5)l-ZP(8,J)>l'(X(2,J)*X(l,J)-(X(6,J)/1.25)*( 1.0-X{l,d) 
1'X(2,J)-X(4,J)))-ZP{3,J3*X{4,J)»XC1,J))+HP(J) ^ ^ 

fiPCJ) = HP(J)+V(l,J)4( (4.0=CZ(3,J+U*C! .0-XCl ,J)-X(2,J3-X(4,a)) 
l*GT/0,5l-2,0*X(i,J3*Z(S,J + l)*GT-l-2.0*Xt2 
l,J)*Z(7,J+l)*GT-(Z(7,J+l)*XCl,J)^WT/2,5) 

U CZ(8,J + U + (1.0-X(1 ,JJ-X(2,J)-X(4,J))^4T/1,25)) 

HP{J) = HP(J)+V(2, J)»(ZPt4,J3*XCl .J)4-2.0*ZP(5,J)4X(l,J)*Xt5-J) + 
lZP(b,j)*XU ,J)^X(1 ,J)-7,p{7,O3=*=(2.0*X(2,J)*X(5,J)-(X(l jfx(b,J 
n/2.5) )-ZP(8,J3=i=(X(2,J3^XCl,J3-CX(6,J)/1.25}»(l,0-X(2,v3)-X(l,J3 
1-X(4,J) 3 )-fZP(9, J3^(EGD-X(1 , J)-X(2, J3-X{4,a3 3/2,03 
riP(J) = HP(J3+V(2,J3*(2.0»X(l,J3*Z(5, J+1 )*GT-2.0=rX(2,J3* 

IZ(7,J+J 3*GT+(2(7,J+13»X(! ,J3*WT/2,53+(Z(8,a+13* 
l(l,0-X(l,a)-X(2,J3-XC4.J3 3/1.25)*WT3, ^ . 

HP(53=HP(J3+V(3,J)»(2,0»ZPC5,J3*X(l,j3*XC5,J3+ZP(6,J3»Xtl ,05 
1 *X(1 ,J3 3 

hP(J)=HP(J) + V(3,J3*C2,0*X(l.J3*Z(5,J+l)*GT3 , 

HP( J)=HP( J3 +V(4, J3*C (ZP(9, J3/2.0) ♦(EGO-X( 1 ,J3-XC2,J3-XC4.J1 3 
l-ZP(3,J3*X(4,J3»X(l ,J)) 

CONTINUE 

RETURN 

^ f » t ♦ ;<i ♦ 4; t + 

THIS SUPHOUTIME INTEGRATES THE INTEGRAL TERM IN THE DBjKrTI- 
VE FLINCTIOM USING SIMPSONvS RULE 
subroutine SIMPSN' 

UIMEUSION X(7, 20013 ,SX(6, 20013 

CnMMON/Al/SSX(b3/A2/X 

C0MM0N/A9/EGI) 

DO 1008 K=l,4 

SSX(K3=0.0 

continue 

DO 666 1=2, 4, 2 
UO 99 J=3,iR99,2 

SX(I,J3=SSX(I3+2.0*(XCI,J3*XCI,J33 

SSX(I)=SX(I,J) 

continue 

sx{iia3=ii^ci)+^?9»(x(i,J3*x(i,J33 

SSXCI)=SX(1,J3 



onof j 


102 


6 6 fo 

!i000 

5uu 1 


C — 


701 


702 


707 


704 


COfiTINUF. 

CnwTlMilF, 

on sono 0=3,1099,2 

n j>J)=bSX(l)+2,0=FCU;on/X(l ,J))-15,0)>l'{CEG0/X(l,J))-lS.o) 

O' X C, 4 J ■* S X C I, f iJ j 

COwTlNUli 

DO 5001 J=2,20u0,2 

SX(1, J3=5SX(l) + 4.0*((EGn/X(l,J))-15,0)*(CEGO/XCl ,vT))-15.0) 
>SSXCl )=.SX(l,ij) 

COriTINUE 

S.5X(l) = (1.0/b003.0)*(SSXCl) + cr (EGn/XCI,l))-15,0)**21 + 

1 (((EGU/XC1,2001 J)-lb.y3*t2) 1 

HFTHOO 

Efll) 

1 1^4 ^*4* ****** ****^**^*^$*iHf*t$*******’if'*^**t***** 
THIS SUBHOUTIME CRtiCULATpS THE aMDUT^T OF VQLATIf.E COMP- 
UHEOTSCFTHXLENE glycol, Water, ACET aLDEHYDElFGASHED FROM 
THE REACITON MASS 
SUBROUTINE FLASH(DIST,NOX,T, J) 

DIMEMSICm X0RG(7),XIC7),F(?) ,AK(7),FnEW(7) 

OIMRNSION VW(2001) ,VDMT C 2001 J,VG( 2001) ,VT( 2001 ) 


C(!HM0N/P1/XI,X0RG 

common'/fi/ak.f 

CnflMUN/FFl/PT,!! 

C0MM0n/a4/VW,V0Mf ,\/G,VI 

Tf OT^XORGCl 7+ XnRGTJT+MRGCT) + XORG ( 4 ) )'/ 27o 
DO 701 1=5,7 
FTuT=FTOT+XnRG(I) 
continue 

DO 702 1=1,7 
F(I)=XORG(I)/FTOT 
AKtI)=0,0 
CONTINUE 



1(3.05143E-17)+(T*»6.0) )} 

PT1 = C (F(5)<'AKC5) +F(6)*AK(6)+FC7)*aiI(7) )) 
IFCPTl.LT.PTJ GOTO 703 
CALI. CUN V( ROOT ,NOX,J) 

IFCUUX.GT, 10) GOTO 710 
DO 704 1=1,7 

FNEW(T3=F(n/(1.0+H00T»((AKCI)»EXPCl .0 + HJ/PT)-l,O)) 

continue 


705 

119 


70 3 


F.NEWT=tl ,0-PUOT)*FTOT 

DTST=FTOT*POOT 

00 705 1=1,7 

XI(I)=(FNEW(I)»FNEWT) 

IF(XICI).GE.O.O) GOTO 705 

continue 

FORMAT (//) 

N0X=1 

IFC CNOX-MOX) ,EQ,0) GOTO 706 


M=y 

GOTO 706 
NOX=0 
GOTO 706 
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7 (/ f! 

7 1(* 
711 

IVli) 

/t »9 


'** VAC. AFPLIKD AT J= * , I*; , ' » * * ) 


FriRWATH A, ' 

»i>u»<;-,/, -CHECK) J=-,I4) 

RFTt'P ''3 ■ - 

tCtiD 

In USING NEWTONS METHOD 
SUhRUU’ITNE CONVCFnnT,NOX,J) 

DIMEiviSiON AKC7J,F(7) 


Ca,'lMUN/Fl/AK,F 

C'HM(iN/FH/PT,H 




35 


1 too 


oO 


c 


H = 0 ,O 0001 ' — — — '■"' — — ™- 

[TA = U 
FflAR = O.U 
ljFlJNP = n,0 
nu tiou 1 = 5,7 

f'”f^U=FUNR-KAKCI)*F(X) )/(PT*(l .0 + P.*({AK(I)/pT)-l,0))) 
Cn'iTINTlE 

funr=funr-i .0 

IFCAb5CFUNR).LE,i.0E-06) GO TO 80 

IF(ITN.GT.20) Gfl TU 60 

R=f..-FUNP/DFUNR 

1 TM = 1 TM-H 

GO TU 35 

TYPE 40 

FDHMAT(15X, 'ITEIRATIUN DOES NOT CONVFRGR'l 

IF(P.GT.O,y)ROnT=R 

CUlviTINUE 


NOXs:l 

RETURN 

END 



2. COMPUTER PROGRAM FOR THE SECOHI) VARIiTIOR METHOD 



jt: li t j)- ^ 'f # I 4 4 

f. 't I 'i* >ii 4 4 t t ^ " 




. 4 f # 


P'rl RfArfflRS 


: ^ ^ 4 t f. f t *1 ^ f I ^ f ♦ # ♦ 


H'ls P.ii'GPA?-’ THfi SFCHa') VAPIAfTU?: •'t.T.if'! 

n'-'lT-fi.!, TF'-lPtRATljR}*: PFOF I L&S T ' h 'jA'if’ < p 

1 :!S ' sTiifiL) JS TAPhIKD n,j flL! TdF nhjrjrt T /F 
'.mT CnT^UC]] FfitlM TTFRAri!,.?- T.J 


t' T ..T' i ■' >'. 

n * 


re :?rArjf’ 
r-M' c^jr- . 


J''..'-' fip'*’! 'AT TL.Hp>f^,<ATnKP t^KflFT.^r ■‘MTATiFi/ Fi-i '»' r'lfiM-r 
IK n^^'A’iT,i!.’ 'FJ’ini; l% /^C ._>• r..T-ilAf, Pj-i.*' 5^. 

I'liiA ■‘'^rTunf) 


.Hi .• i: } 4 » t 


i' 1. 1 4 


P^uGRA.1 


f**4 + t* 




f. 

r 

I 

> 

ic. 

I 

¥ 

u- 

f 

$ 

f 

i 

:|c 

♦ 

♦ 

t 

% 

t 

I 

f 

♦ 


X 

4ri 

4 


f\Cu 

I>T 

'T 

I:;--'- 

K 

ZP 

Vb:,T 

(/G 

V T 

Z1 

PI 

VP 

EPoT. 


EOpT 

ak 


DPBAR 

P 

y 

OEELX 


STATl- VARTABbES 
aDJOIET V/^HI ARuKo 
ACTTVATTUM E iFKG’i 
l-REnnEfiry t-ACTOP 

RATE rOPSTAETS FOR ThF pMCTIO-A T O'UFyF 
CnwriFOF.BTT'ib STAbF OF PFj FOk-ATIO.. 

'•'TOTH OF THP Tf'V, IF'iFHVAf., 05FiJ V T 


I’ 




TOTAL TI-iF rjF i tVlA, 


'jOKA !’T0'. .'3 
IF EG GPLi'4' f ! TrrF r- F ^ 


UJITIAL COuCE'i rRATirp 
mTAI. TRESoUHE, AT'1. 

RFArtPR 'lE IPERATijRE.uFviRFr-' «-'Ei.. T, 

’inoPERaTORF At 10 .1 PuT*. rj ,r TIME 
TFMPEPATli’’fc AT ?uOl PPIOl.'; np r!.,''‘ 

U'''IVi:P;,Ab GA,S COORTAaT 
PERT V ATT JE OF L -v.H.T, 

MOLAR VOLLME OF UMT 
.-lOliAR vnL'PP'' OF GT.TroT 
''•OLAP VOJjHhF of .^ATEk 

TnXAL V()L'JmF of the RFACI’TO': ^ AFG , bl’i PbS 
Z OTVIOFiJ AY VT 
oRJFCnVfc: KOoCTTON 
LERTVAt'TVF: OF VT W.P.T, f 

EPSILO.^-FAC TjP PETfRaIPITG ThF ««G.an.i|.>P ot 
TAKFR JW T.iF rijPRFCTlOu .OF 'r’FOn’FkAIPKF *='-<0--'. 

043 ITbPATUO 
nPTTf.iSM LPSTLO-^ 

VApnpF rRrsSijPE , atm 

AI.PHAl - MEIGHribG PaPAMF.TEF FQP TiiF riFSlEtD Cu 
ATint* OF DIETHYLF4E GLYCOL TP X‘4: POLY.IEF 
ALPHA2 - 'VEIGHTTPG PAPA'^ETSP FUR THF PESIREP Yh*. <•!• 
UPfiAU OF PilF PlJLY**EP 

ALPHA-1 - WETGHXImG PARA'IETSR FUR t’HF FUPAATI04 ‘-'h 
3I0E PRODUCTS IN THE PEAC’^a^’ MASS 

IkUMRER-AVEKAGF MOT^ECULAR height of Trip POT.YPER ,, 
VARTAHLF TOVOLVKO £M TXE RirCATt ro?lATin.4, THIS IS 
A { 4X4 ) 0 A 1 RtX 

variable TaVULVFD IN TOF 0 RyUATlfiiS.THTo IF A («X*> 
VECTOR function OF TIME. , . ; . 

PEKTUR&ATTUM TN THE STATE VAHJaHLFS 


STcF ■■ T" iiF; 

I i 4 ' 4 p r.V .1 *» 



1 On • 


f i J. « » ^ 

(t 


■1T 

••r 

I- M * t i 


.•'T.-'f=r OOi-'TV.VrtVF nf- i'Ht hA ...H.T, 

.,.FC'',,n nE'U’/4T[V£ .)F H A'' I! £<’ • T 

I”' I VKRSAij 3AS Cn.45!rA ;'»’ 


T ;£■ ' 1 = 
. iF 


S; t * ♦ ^ ^ t t ^ 

r<~! 
C*'> 
CO. 
C 0 * 
Cf- 


7,2Onji,v(U^0ui),^p(a,iun 
t'f J , f'JE J(lnt ) .IPduU , O ( 7 '» tn 
/.f -1.2 0)1 ),/,iri,2uoi),7:i( n,EHf.n 

JOI) ,/0-' l'(20t I ),v,'}C2<-l),7t,(^...jj j 

i>ph A H { 2 0 n ) , EP c 1 s 1 , t ( 1 s j , }■: 1 ( I ‘n 
Of J, Ul01),nKTp(l,4,iOl) 

.)Pi.T,x f i , 1 , 1 'll ■» , OPT Cl , d , f i , 4 , 1 '■ n , '‘O c 1 , * ) 
Mr2non,Px(4,n,XE( 4,4) 

*1**4 t 

"tr /.M /.SSX(t>)/A 2 /X/A J /7 

.’•.•i'x/A 4/t/'M,v!jMr,7c, 7r/fc:i/Eu,an,Pv; 

•'■I* /hb/ 0 ^ /An/T'J/AV/V 
op-/AP/i,r,u 
i;r'/FFl/PT, 4 
!iJ ./Pa/i)JC>fr' 


-ll 

t . I ■ 

.6 1’*. 

: r. 

'to ' 

.S 10. 



'61" 

1 T 

. t”:’ 

SI''- 


1 ‘ 

■S (f. 

■ ) ! 



uT 

! 1 ’ 

-SII’ 

>■*1 

s !!'’ 

-6 in*. 

•'U 

■ 1 0* 

..SIC. 

l’» f 

-fi;: 

. S 1 1 > 


CO ■•■iiii'/Ai2/r''pr- AP 
cn iMh'./P 11/p/p 15/nFTp 
CO !MU ./AIO/P/A 10/012 
CO : ^.1 )/A19/FX, aF 
CO ''■•j';/A24/0/A2'i/nCLbX 

■ 4 4*44 4t****n 4* 4 444*4444 -f************ 4*4** * 4 ♦ * 4 » 

ORAL MO 
EC3=1,'( 
ii=-1.0 ^ 

0IMTi:C40,540) 

TYPE 540 

540 F0RMAT(1 X, 'Cavp; I’HS VAbOES Of A0P0A1 , AuPHA/ , AOPFA? T.i 

IDOUBLE LTOFS' ) 

REAi)C5,*3 .ALPHA! , f LiPHA2,ALFHA3 
.'(RITE(40,33 , ALPHA! jALPHA2,AI.POa 3 
: i PORMAT(U,'ALPdAl= ,P!5,5,2X,'AT.PHA?=',F1‘^.5 ,')x, 

1F15.5) 

! H=101 

WB!TE(40,673) 

; TYPE 673 

i 67 3 FOH.maTC! X, '01 ve rOTAL pRESS'lKf') 

! HEA0C5,0')£P'r 

! 4niTCC40,573),PT 

1573 F0HMATcU,'PT=',F15.S) 

I UPE‘'KUfiTTs23) 

I HrAO(23,54y) . (THEWC n ,J=) ,H) 

I 540 FCIPHIVTC il InlF? jfi) 

* 1213 rn.PMAThxPcUES.S TE^p pPijFILE FOk 1’PE SFCn.^,P ORDEF 
j imiTE(40, 12437 

WB TTEC 40,549) , CTWE^j(J> ,J=t , 101 ) 

6=1,0/1000.0 

RG=l,987 

H0=0 ' , 

' E0{3) = l8500,00;FO{4) = 298u0,0t) 

; Eol5)=298Oy:O0;EOC6)=29800.O0 

i efn7)»l7609,00;FO(8)=17600.00 


1 f # 


♦ t i 




f 

# 


* ^ * S ' ) 



f n ) 


t.. ■( >:■>. ;n 

;;■ ( }) = n.i3;i,n(4) = j7.«.4ixnisj = i7.«; ^ 

;v = ( n J = 1 7 . ^ 4 ? Z n ( 7 > = 13 . >5 5 5 ; 7 0 ( H ) = i 3 . b 
?/’C‘Ti = ?2.044 

.Ji !=:‘.j.J+ I 

1 TciC -tti , t ,ZC; 
r7.^>: Ki.'vH 
C •’•v'-'.'-.l ( U , 'Mu=' , 17) 

.'■fa c-MMh) 

..P ITi' c 4U , 1 ♦>) 
h;!.-- 3TC1 X,l?0nH„)) 

'ini'? PAPt CALCjT.ATi;:r> T.lF '■'uLAP VJl.i 7*.;.^ A in p, 
M r POir.'TF OF r£:-lF 

i>n iu i=j,9 

4( i , t Is^i'.O'frLXPC-F.nc n/cnrt'T/.Fi^t i n+zn( I j ■» 
ZPU,1 J=feO£T j*z( I , 1 )/(Kn»T n ) 

C'l i.’IMSt- 

r^(u=T^zwu) 

V 0 ’ ! '«' ( 1) = 1 <) 1 . 5 ♦ c 1 . 0 + 0 , 0 i* 1 4 ♦ ( T C n - 4 I 3 , < • ) ) / U' fVi 

\?G( 1 )=rjO,6+C 1 .p' + u,001 4*{ IMC 1 )-a1 1 

i/-Ul ) = Cl9.472 + «,'.,i26’t'(T7f n“413,r))/i0fii ,0 

i I 1 *** ' ' 

ixi'l 4 d = J ,100 
U=U-) 

Dfi 1 ♦ 1 = 1 , 20 

f f n = ( i v&iric j+t )-rF£A'i j) '»»D/2y,o+T’'!ew( jj 
1)1' 15 i 1 = 3.9 

Z( It ,J1 )=#iO.-}*Kypf-F.|f It l/CRy*! (1) ) + /,Dt r n J 

Cf'Pt'IfJUF 

vr(JI)=t(t) 

v(>,/r(ji)s:i9t ,5*f t. o+u.poMtCTM jn-ni.ojt/i 
VG( Jt )=60,b*(. 1 .0 + 0.f’014*t T,a J1 ) - it 3,0 ) ) /lO K 
VW(.lt ) = f 1 0,4 72 + n,o^5y(TM(J1 )-4l 3,X ) )/1uftv>.o 

ii = l)+l ,0 

ji = n + i 

lOll'UL Cr‘4)lTinH.n 0,M int state 7AKTA*^uFb 
X(l,l)=l .0 

x(2,i)=y.o 

A(3,l )=0.0 
x(4,i)=0,0 

X(5,l)=0.0 
XC6,l)=o.U 
xc7,n=o,o 
CALL RIJ.^GI (A) 

FnRFATC?X, 'X is: ' ,2tet 0.3) 

WRITE { 40 , 6 7 R ) , ( ( X ( L , J ) ,T=:1 , 2i>(a , I Oo ) , 1 = t , 7 j 

WRITF.(40, 1 5) 

fYpF 791 .'■"JfI'TSTK 

FORMAT! 1 X. 'FUR 00= ' , 1 2 , ' AMOOfT ^L ASHEOs ' , F 1 5 

«RlTE(4y .791 ) ,!'in,niSrF 

OlSTf=0.5 

P'INAL CO.vDltlOMS 0-M THE AOJfUNI VftftlABLB:5 
VCi,2CCl )=n.o 

V('2,20Cn = 0.v 

7(3,2001 ) = A1.PF At 
V(4,20<'1 ) = 0.0 

14P=( X(l. 2001 ) + Xf 2,200 t)4-X(3,200n+K(4,200!)> 
VT(20(;n = Vt)A'T{200f ) yMP + VGC 200 1 ) *X( 3, 20ol J + V W 
Z1 C4,2('0n=7t4,2f0l) 

21(9, 2001)=Zl9, 2001) 


LlTlOj 


state '/AkTARuES 


'AMOOfT plASHEP=' ,F15.hj 
AOJfUNI VmihBhKS 


00! )>/2,0 

1 J + VW(?U(<1 )'»'X(R,XpO J > 



lOd 


^ f't. c n 9) ),',rrfn i m 
■/.1 f f ,7 )=:(Z{ l,2CUl )/VT(?OOn ) 

• < ' C - ' . ” [ - UiS ' ^ . . 

:«■' f -J-t, 1 t. V ■ . 

->' i rbf 4''' ,h,)0) , c ( / f I ,j ) ,j=i ,2001 , icb) , rsi ,4) 

■i < r>f I 2 # ..0’ 

1 .d: ' , 21 >%n .3) . 

Ci'l.i, 

rts^bPi.M »x( 3,2O0t ) + ( AijOHA2/2?3,0 ) ♦';sy 1 1 j+AbpM .3^r:;.':xf 3 ) 

I + :> S 3 { 1 ) j 
■ lYiO ': 2 b , fl 

I Oi ■■\-\ ( J X, 'PJ = ' ,F7'‘.7 ) 

. 0 1 n .CAl. ,lr>l 

*.'p ( I’K(40,23) ,P1 

IVI-K 113 

,C1 

If'C''! .F(!, ' r JfbiT.l 4? 

ITv .S jCuTil i? 

C i-’T i/,ij COhrUTirN,'! nj TIIK "P" i)TFFFP£Hr?A], PyUATinr.S 

!.>!• 70 J = l, t 
liit yu (<. = 1,4 
Ff 1,1 rlODsO.O 

'r* rof Ti^HF • 

i'l'jAh COwn iTinrjs THt; "«>” DifPfc:;»(-:*’rTA?. Fv"AT,i.nbiS 

un 79 K = i , ^ 

i=j 

n ( o -., 1 , 1 ( J 1 ) = 0.0 
79 CfJ'JT 

C f‘ir»'iAi, cn.jniTifi,,F I’Mh pfrt.jphai ru‘! PvUA'^ini.s nr i'-'r. 

C ^STATF UAPTAPlFS 

tjfl P'j nri , 4 

• 1=1 

I)CLLX{K,I,1)=0.0 
ilp Cfit' TIf'UF 

CALL ALPHAS, At, FHA3,AJ 

VYpr. 1600, C((P{T,K,J|,J = l, 101,101, 1 = 1, 4) 

*/niTF,(40, 161 


13</00 

1600 
17 00 . 
IBOO 
13001 


2929 


CALI, RUFG4 (ALPHA2,ALPHA3, A) 

TYPE 1700, t (0(1,1 ,0) ,J = 1 ,101,11.) ,1 = 1 ,4) 

TYPe 113 

RFAi >( 5 , 45 A ) -hi 

iF(FJ noou 

GO in 13001 

CAi.l. HUAGSl alphas, A f.PHAl, A) 

TYPF lHC(',( (!)FLI,Xf I ,l,vJ>,J = 1 ,101,10),T = 1 ,4) 
FnP 0 AT ( U,'P TS :', lleU , i ) 

FGP'IA? (1 X , *0 TS: ' , I IK! 1 .3) 

FG'OGiTUX, 'O fILX 1 S : ' , 1 1 K 1 1 . 3 ) 


FG'OGiTUX, 'O fILX 1 S : ' , 1 1 K 1 1 . 3 ) 

TOIS PART CALCUI-ATKS THF OPTIMUM FPSlT.j'^ 

TYPK 26?9 

O') 601 K2 = l,lo 

EP(K2)=0.0 

ai ( K 2 )= c.P 

Cohtioof ■■■:■ 

l'lipMAT(ix,'c;ivfc: nie yaguf of si’AfttiNG epstlon') 



fc 1 5.5 1 


■ U > )-''r‘' 5 u 

■fi-'Af t,7 x-iji ,rifr;;i) 
f ' ■ ■'. * 1 1 A, , ' 1 j)K p’l r ht j Ni, tPs Ti.rjs ' , 

. ‘ i ?•*.( in ,7 .u-, j) ,.n j 

5151 = 

,J = .J 1 

^ = l 

.>'» 52 1=1 , 1 

I ,*■ J=-,>f l,h ,J) 

■ , lj=! (K, I tvO 

n=nrTi'(H , I ,j5 

rt’-ijf. 

I'MJ'inijTC.’V I’Tfjf f'r TH!-; .l^ VP'lC.Ff^ " >R'‘ /VJ-) "Di-'T" 

„r'rT = 0 ,;; . : ' ' ■■ ■; ■■ 

O' I ?! 1 = 1,4 

n-,)!* rs ivif- I + 1»F1 (f , n *0!Jf 1 ,rvi 

Cf' < 1' I Mir 

iMivtTciLlCAfTu' OF t'f!F nP I’PICFS ”KP'' AWu ‘'Dn.LuX 

? I ) »•• ( I li >, s 5 , 0 

?■” ?2 i = l , t 

Hr!K!,bX=Pi!r:blX4i:SHK ,T)^r)FLf,x«U ,K 1 
CC'rTl'.tlF' 

rFK'M = CKTt uU+DrFT)«-EPSr.’FHT?tJl 
TFK'.«? = .l‘f-/{ JlxnDKi.! X 
i;i((il)=Ti4F,F( J j-TFH*!J -TFW''5 2 
rr«(,i)r + n-K‘U +TFRV 2 

rF(TH(,n .r.T.poo.oKiofn 5 
rf'CJ)=oOfi.o 

lF(TK(J).GT.40o.U)G0T0 5l5l 
XP{J) = -10(..0 
CTiNTimiF. 
t>0 60 1 = 3,9 

'6 C 1 , i ) = 6 0 . 0 * F> F ( -F. i ( I ) / ( RG’K Ti? O ) ) + iO ( T j ) 
CnfiTIwijF 

V P f T I u = 1 9 1 , 5 ♦ (1 , 0 +0 . 00 1 4 ♦ ( TR ( 1 ) -4 I ;u 0 J ) / 1 v 0- 1 . 
VG ( J ) = 60 , 6 =F ( 1 . 04 r. . 00 1 4» C TR (1 ) -4 1 3 , t) ) J / 1 0 io . 0 
V2(l ) = tl 9.4?2 40,U25»(T!<(l 1-413.<'J )/1 0»>0,«;. 

Jl=2 

TOID-TPII) 

DU 152 J=1 ,100 
1>=1.0 

U!! 152 1 = 1 ,?(’■ 

TCx) = (xfU J4l i-rr c.T) >ii'n/?0.04TP(J) 

D(j 5 3 11 = 3,0 

zcn ,wn)=60,0*EXPf-EU(il )/{RG*=rCI))+Zn(Tl)) 

CUf'TiWUE 

T.'J(JU=T(1) 

VDhT(JU = lMl ,5»c 1.0+0,001 1-41 3,«1) 1/1 li'l 
\/G(Jll=60.o*h .0 + 0.001 45 J -4 1 i , 0 ) J / 1 OoO . 0 

VU( J1 J = { 19, 427 + ‘i,025 + (TF(jn-41 3.0)5 /liiOU.O 
U=D+1 ,0 
Ji=vn + i 
CUHtinuf 

CALI. RUFJGl (A) . 





llo 



r' n 1 “ ' 

\ '’(. 1 1 , 'Oil, j r-sr 


•r-ypr 

'^7^5, f (.XM 

$s 

r-\i J. 

OJO'l’O--! 

tf' 

f , r 

= >0 1 iinl ♦y{ 1,2 

ll 


( '1 ) J 

I ' 

r Y p i: 

1 '/.i^'O 

, / 

'.Pi r 
1, 

r’fif '■ 

1 )=:«Oo 5 l 

B .! 

■\”t 12, 

;•> 


(5 ,4 LS j ,h 1 

i 

1 F(*^ 

i.l P. '.s'iCOTu 


j j ( r 

1 ) = PT f‘FW 

, 

* 

r/i'F 

7 1 ,[.PIK5. ) ,0H 


, i> j f 1 ) + in i t. t 


' i .-). «; 4- * 


72 


!V!>r H/ 

,i'.l 

ir(-'i.hU. * 1 .'' jr,nT’;i 72 

KPSl, = k;pPL/2.U 
k 1 =ri ♦ 1 

IF IM .i.n'.l .UGnrn 57 


: ^ t * * f- * r 


npn’.runpf //,. ) :' ) 


7v typf. ) 2. no 

TVPF 1 ; ^ 

125 1 X, »p[) xoli ,-.i\;iT Tij GO BACK IB 

-a- A 0 t 5 , 4 b 5 i .Hi 

i:U 

m! ITGCi.O, t u) 

‘ 1:7 01 , i?,xn 
GRTTt.( 40 , Itil 
«!aTE;f 40 ,fh) 
wRITc,(oO, 6 b) 

TVi't b 5 

o5 Kr'Ftf.'ATd&X, '1',3 hX, 'E'} 

I"' 7 ': t' 1 = 1 , 1 

.HUTl { tu, 7 i ) ,01 ( r 1 > .EBIM > 
wwrre( 60 , 7 n, 0 Hi‘i ) .EPtKi) 

71 FnH!<AT( 7 )i ,E 17 . 7 , 3 X,»^ 17 , 7 I 

7<! Cf'CTI'.OJt: ■ ^ . , ' ■ 

64 FnRHAT(,?X, 'FUPT=' ,K 70 . 10 ) 

TYPE K.'i'b 

1008 F 0 ( '-ATtU , 'Gjvt; XHE VALUE UF FOPf') 

RF-:AP( 5 ,’*) .AT 

fOJPTsAt,. ' . 

.IGl 

WRIT Ef 40 , 1 0 ) 
v/B 1 TE( 40 >, 64 ) .fitn 

*C THTS PAF'i rrpPETTS the old T'EBPFHATnRF PHnpiLF U6Tr,r, 

C THE (iPlIEHn F.PSTLON 

DO 5152 tJl = 1 .101 
.)=J 1 
K=1 

DO 5252 1 = 1,4 
0 R(l.K)= 0 f I .K.J) 

KP{F,T 3 =F((< , 1 /J) 

DFl,LXRtI,E)=DELDX(l,K.D) 

0 FT(K.I)=DFTP(K, 1 ,J) 
i'525? T 1 ’*HiF 

PC HHDTXPDICATTDN OF THE MATRICES "OR** AFD "OFT" 

DOFT= 0,0 

DO 2121 1 = 1,4 


nic 


ly.i 




+ tF'j ) 

THH ^.F rRircp "kr" 


»f.u .’•rcFivXh” 


; ? 
r -.1 

. ' s- 1 -11 




; 1 ? 0 


115H 


455> 

lU 

42 


»•»->■* t-’* * 


46 

45 


I’//'.* 

’■r.''" iriTCATTi-.f! np 
!’»i ' y y. i =: 1 4 

r.Lx=n.Ff;f a'+k?h»^ , r i^ufi.lxrc i , *■' j 

f! i ' 1 T i i 1 1 1 F’ 

■i’n 'n ! (.1)'+nFFT) ♦KnFT*MT2(.J) 
i'f, .'•■'k'shi Ff J)-H;i)ruLA 

; * { v!)s t f.t !- (il) - I KH»' l-TrK-’2 

( J4= I M y + +TFF*-’V 

I K( •! ti( jT .r i rn i 

'■ r (.1 )=i>F( , 'i 

JFCi’HC J4 .GT,4no.0')GnTn 5152 

t:(.i T ! i-iiF’ 

S 1 1. 

f)I’ 02U G=l,lb1 

CFi.TlMa- 
wF 1 '!>. f 4 0, 1 PI 
* F ITLtoC, 1 6) 

TYt'K 1 2F,P1 ,{■>!('{’. 

FIlF-'iATCl X, 'P1 = %K1 5.H ,5X, ,F1,5,«) 

TYFF i 13 
kt’AiH5,4f„5) ,H1 
U‘tn,Ff>,»y'jGtni) iis9 

v'.riT'f 4 2 

ir{(4B5(i'T-Fi.'h:F)).r.s.o.ooui jCATu «? . 

IF( JJ.GT.S/GfiTU 42 

fYFfc: 4b, (TwFwf j) ,vi=l ,ioi j 

vKlTbl 4(»,45) , (TFKWCJ) ,J=1 ,t-l) 

■BMTKC li , 1 o) 

FOpy ,G 1 I A n 

FOHMATdx, *ri> yni( to PRf»rbFu( y/.-.H ' i 

GO'IO 11 

TYPl 4b, tT,\F>(J) ,J=1.1 pU 
.'fUTb(pO,4 5),Cff3E^CJ).J = l.lOn 
‘'iRITh,(52,S4ti) , tTKF'^(0 5 ,J=l ,5^1) 
on 46 J=1 , 10 I 
lF(f.jF/j(05 .lT.600.t>)G0T0 4 
TtlFiiCJj5:t)00.0 

IFlTot v;(Jl .Cl .400.C)G0Tn 46 
fbK--l( J Js-iOO.O 
cnriTi FuF 

Fnp'/.ATClX, 'Ti-Kv; 15 ; ' , 1 01 F7 . 2 ) 

.vHl i'h’tAO , 1 h) . , , . 

OP(jAK(2''5i J = (FuO/(X(l ,2001 )+X(2,2Ui;1 )+a( 3,2'^01 J + X(4,20u 3 >) ‘ 
pPTTFCbO ,P90) , (Xtl ,J) ,J = A ,2001 ,2ul 
.'iRlTEf 50 , 1 6) 

WRlTfc;(5';;, S49) , (X£?,.T1 ,J = l ,2001 ,20V 

WRlTlfso'HOO) , CXt3,J) ,J=1 ,2001 ,20) 

HRlTK(50,lb) 

wnTTb:(bO,«99) , (XC4,J) ,J = 1 ,2001 ,20) 

«Rirt;C5o, lo) . , « ^ 

WRITECbO, 599) , (XC5,J) ,J=1 ,2001 ,?U) 

WRITL’CSO, 16) ^ , 

wrUTb:(50,By9) , CX(6, J) ,jsi ,200! ,20) 

WRITE £50,16) 


WRI1*EC50,899) , tX£7, J) ,J=1 ,2001 ,20) 


112 


• ■ f ; 1 o) 

'■'ll » If^l^RnPCI) , [ = 1 ,2fMn ,?.>■» 

■^•1 t ',f 1 c, ) 

■•.•iT< 1 A, 'X T^-,: ',iyiFl5.6j 
* 'I rt 1 > , 'rrr.Hf-- ts; ' , i oiF]l 5,6 j 
i' ' •’ I ! ( . ( t , 1 2 ) , h ! I 

.STTiT’ ■ ■ 


12 

4-444 


IIU 


2222: 


J , v^'l 

, f 2 


T”T'^"”sIi~7nTjTT7r; T”TR77A'T'Ks'’TTtr”5 

(<■>1 f OliFJ'H-nROFP PlJfiGF-KUTi A >‘1 

lit I'T -.F PUfH.J lAI 

-‘li .; iSlfs X{7,2n05 },xnP,'U7i,Cf l7i,CPif7),rP2(7j, t?.? 
(f rr3( 71 ,rp6(7) Jt ,'nf ^,^..1 J ,,r^ 

. F..‘.-.srf'! XP(7) 1 ) ,V||- Ti7i;ni ) , - f 2 

iT FnSlP., XP^Fj (7 ) ,7| (y) ,F(i{9) ,ZuC71 
CP 'Pi J/A;'/X/i,4/Vw,vn.lT,VG, vr/R 
pf) ' 1 i; /p(,^;l;1 i/Fp,;{);/ci/'Kij,znrr'G 
pn^.'ii’Pi/Pl / X T , XfiT'G ■ 

CO' 

CP..''*P‘'/r[- 1 /PT,y 
CfX'DufjyAip/in-Pnr 

;< F HU P ■ 

PFAG ''PP 
A=A/!i.P 

;>'■ i!ip .■n = i,2P<-F 
j=ji 

.iP=(A(l fVi')*Vl2,,1J + X(i,J)tXt4,J))/7.<-| 

!lPi,AH(.n = ( iFcn/Ppl/'^.O) 

V'’T’(.J jsvPMTt j)#MP + VG(j)*X(S,,1) + vWCJ)*Xfb> J) 
pn 99 7=1 ,7 
Xni^GI f i)=X(N,,T,l 

ct'r; riuUE 
9=1 
P=U , ( 

un 1 2 1=3,9 
zja )=z(i ,j) 

rj\ ii'i'-JUF 
C!) rriPuF 
on iol 9=1,7 
XMNG(r))=xru,J) 

XPC' j=X(.<,d) 

CnnTIPUf. 

T1:=T'i{J)+P'««CT''l(.1 + n-T1iC.1) ) 

Xl=191.5'«'fl.U + 0.C914*CTl-4n.(,-j)/tU0t.*.F 

X2 = 60.p1( 1 .fMU,O014’»‘CTl-413,(i) j/lOOA.P 

X3=:( 1 9,42 2 + 0. 925 1( n-41 3.0) l/loOn.C 
MPP=(Xi.lWG(n+Xa9Gf2)+Xu9G(iUX0»G(n)/2,'.' 
t/TO=Xl ♦■Fi’P + X2*Xn9'U5) + X3*XURG(6) 
i»n 1! 1111 = 3, 9 

X( I ,J)=00.01tXPF-Fuf I)/(HG*Tn+ZU( D) 

CnuTlOUF 
on 22222 1=3,9 

IF( (i.F0,4),i.R.f 1.ii;0.9)lGnTn 22227 
XI ( r,J)=zf 1,9 wrrt 
CO JTIMIIF 
2.U4,JJ=Z( 4,.J) 

ZU9f02=ZC9,03 

U=li + 0.2 

CALL EyTNKJ) 


' ' I .* . > J f 
7 . ■ I 11 

I'l 1 • I I 



O'- 1 ( 1 ) = i 1 

0 ’-1,1, UD 

.'■) S-'B 1 = 1,7 

i)=:A-frp( n 

Ai (rjs,.iiM'r(T) + 0.h*CP2(n 

ili'-. 'Vi •)-]!<' 

).'■ 1)f, 1 = 1,7 

op jr c )=AirHf n 
V<(T) = xr'Hn(T)+C’^i(I) 

ct'or 1 >(!!•. 

(.1) 

.ri !'>*■) 'T= 1,7 
( 0 = Al‘rPQ) 

“m i? ^ Pr 1 n 'T ^ n +C P 3 ( I ) ) 4 C Po f 1 ) J / O . ') 

X C 1 r *J 1 X t * f X ) -I'P P f L 1 
O.i-TX'MF 
•l = ,' + I 

CuTl) -v44l 
'•’) '447 1 = 1,7 
.v( ! ,.J + 1 l = Xl 1 ,-7j 
X( i = ( T j 

CO '.’f X7!JF 
))0 5 7R .r=B,4 
/.C { , J) = Z1{ L) 

Cni-/rFn!K 
!'=i':(d) 


goto 3 4 -3 4 2 


Cl=ZlCi,J)4XijP(;(l)i‘X0OG( D + fZl ( 7,..J)*xn;>r;( t ji'AnHC(6)/2,5 > 

o=;.t n ,.n ♦ i.o^'i.nt.n-KKRGd )-xof-..f2i-XuP-;f4'> }+ 2 ,oi^zi t , j j i • 

i+2,'^*Zl l,7,J)*KnKCd ) 

(:i = (/a/pj 

C'2 = l:;XP((-td*A1 ) 

X(5,-I + l )=Cli'ti .n-r2)i-X0PG(5)»r2 

cnfitH'.'UF 

Of) 1P7 1 = 1,7 

xriKr,tl)=xf i,j + i) 

C<)7 Pi.iDF 

CAbu FLA;;fUi)i')r,xnx,T,J) 
tF( "iUXJXi.i-) OOTij 100 
<(b,J+l)=Xl (5) 

X{r,,.I+l) = V 1(6) 

XC7,.,i+l)=X,H7) 
i)l S^KsDlSTF +DTST 
C07 rirajF; 

.)=■■) I'S.O 

RtTUH'l 

COLS, 

fHIK S(!PHnijTI.NK Cnf7TAl‘'lS THE SfArE VAPIABLE EniJAiTU’'3 
;>n}5P(iUTT'‘tF; roToi c j) 

OTflKHfiino FP(7),XPC7>,ZUP,20f)l) 

COMMUfi/Ab/'M/PiPiGl l/FP,Xa 

CnoourJ/A9/e.Cij ^ ^ 

FXH1)=-2,0»Z1C i,J)*(XR(l)*XR(U-C(2,0*XR(5))/0 *>>*(1 i l-Xnf 

12)-XIU‘)) ) J-ZU4.-»)*XR(1 ) -2, 0*21 )»XRtr))-2j (6,0 »*XRM i* 

lX!<(l)h:6+Zl{t,J)FC2.0»XP{2)'»'XRC5)-(XKCi)*p{p)/2.6)) 
l-SI t8,J)*(XR(2)»KR(U-CXP(63*(EGO-'XR{l)-XR(2)-XK(4))/l->5 •) 


1 1 4 


l3,vi)*XPt } 

1 )+Zl 


*6^ I » J ) 


n '' ! ! ^ 5 -M 1 ^ 6 5 / 1 , 75 1 1 rh“,.- . p I n -xp ( 3 ) - \* 

!) = 2..)*/J )*X»('i)+Zl (f, jyArtl J 

)=Zi J*{F3n-XR( n-KR(2)-XKf n* 

.j=zio,.U'^('f‘^ri)*XR(i 3 - 4 ,('y xkcsi * f c;f*j-AP( u-xp ( 2 )-a‘ 

.J ^ X ( 2 , 0* XP t 2 ) *XP ( 5 ) - ( X,u n * XP (■ 6 3 /? . ■> ) ) -7 . H b , n « ;>. 

.j=zi r?, j)n:>.uyxHf 2 )*XhC 53 -(XPn )*xw(r )+ 7 j t '.r 

uu-u . 2 .jiycEp.ij-xp(n-xp{ 7 j-xpc i )>) 

') =Z U4, lU *X’’( U+Z Ui, n*XR(l )+ARt 4 J 


i A ’“' I I i 


i >'4 n I I 


A ^ I’ i I I 
i .|, I 1 I m 7,] 
k { I 1 f 

i f i ¥ r, f / 1 


nn.'X MiHf.M)1if)K CX[,C!)l.RTc;S Till- A iR-i'M' uK ’'i -L;' LT c. T 
M'lr’." rSlKl H iT.K'-f-. <;ijK!r.i(,, sAT-^R, ACFrALni-.'n’i'»t,H'LA.<^.,r^ F 

T'h, R7A(M't!i M,\*o 

iHhi’-HITT'fP Pf,ASH(nu?T, i'lK, l.J) 

i) T If'j AS L! i .. :<ivn t 7 J , X T ( 7 J , F ( 7 ) , A K f / ) , F,*Ff, ( 7 ) 

cri'., XM'i/r i / X T , X.'MA’; 

CO !^'i!)’)/Fl/4K,F • 

CO-iMuM/PKI /RTri. 

FTOT=(xnkr,n )+xnK‘.U7)+xnRf:c3) + xnaGt4jl/2.R 
on 7U.t T = S,7 
FTtJTsF rn I + xnuG( T ) 

CO cntnjfj 

on 702 Tsi ,7 

CC t)=Xi)fo;f D/FTcn 

AF(T)=U.U 

cnr;TiMiiF 

■\K(S ) = in, 0 + * C7t ,(,1 - 372 ),C/T- ,.t 12* AT IG I 0 ( ’’• ) V 7r>0 . 0 

, \K{0) = lO,0**(«.00‘U<'i -17 57,8*33 /( 1-3 3.2 /4 J )/7no,.' 
AK(7) = C55.0*tXP( 41 .215 329-(4860.8RS4/ ri-(S.LM!l43*ni, 

i( 3,0*31 4 iF-rnyfr+ys.o))) 

nTl = tCFC5)*AKC5)+FCft)*A'<(6j+F(/3*A)X(7)J3 

IFCPl'l .1/1 .FTJ Gntn 703 

CALO Cn07(KnOTl .kOX.J) 

lF(onX.CT. 10) GOTO 710 

on 704 1=1,7 

F ilFM I T ) =F ( n / { 1 . O + HOuTl » t ( AK( T ) ♦£XF{ 1 . 0 + H )/ Pfl - 1.0) 

CnoTI'liiF 


K'Jt':iT=(1.0-R0nt1 )*FTOT 

l)I,^T = F-fm *HnOTl 

on 705 1=1,7 

XT(I) = (Fi;Fuf ) )=fF!,FwT) 

IFf XI( O .GF.O.O) GOTO ?o5 

COuTIonr 

FOR'-^hTC //•) 

liOX = l 

.IF( (oox-srixi .Fy ,0) Goto 706 

M = il 

GPin 7PG 

!'inx=o 

T () 

FnawATClX, '** VAC. APPbTtr? At’ 

T¥PF 7 1 1 . 

FOHikTt 1 ix , '* SOHKTrtlNG HAvS GONE rtRUMG' , / , 'CHECKI 0= 
IFCOOX.EO.O) GOTO 700 


f 4 t * 
p, 

r>. I ’ 


fi-;c 


f ** ^'1 '"I ■ p ■ 




' n ) 

M f u ,j 1 I ‘I i 


rMflM.-'-inV ft'h KI><ns THt, KfUiT I j-'* ^*VirK> 

rur>v( 

.».r p-’ ,r^io.. Aff(7 ) ,f(71 

f'.V-MtliV*'’}- \/P\ ,M 

■• i )[ I I'sH .t/uO '"' 1 ^ 

LT.,= ,, 

: j F ‘ i i 

i'i liUi) ls:i,,7 

^;5.k f p ■+ p (I ) ) / ( pT ♦ f 1 .C+KHOT* f ( Ar. f I VP i - 1 . > i) J 

aF'jMn=l)!-MVP-»l-AK(.TJ*Ff lV(lA^(J;)/r•i')-^.np/fi T+f ( « 

1 ( V(I j/PT)«i.n))**'>)i 

CPriTlPllF 
F(lMP = Fijf;w- pn 

{r'(ArM‘HF(ifiiO.T,i',,t,0F-0ij) Gum «t 

If- ( r’rr!,f;v, 20 j Gu 'I’u bO 

un{i-l' = P(iU'r-t* nii!V('FU,iP 
lTi=lT,i4i 
ui) Ti) 35 
TV I P 4t) 

FflFMATd GX , ' n'KPATIIlM UfipS *!UT Cu'V '^KGi:: ' ) 
iFiRuO'i:.(,T,u,o>f'>ini:i=ni)Ln' 

C’(P,"'lNi)F 
!UI<=t 
i’rpnpp 

V* ♦ *4. 

'HUS siiHpnHTiijp [s usfu mn bni.vrjr, tmk Arujni'rr v'Ai t 

F.CMiATiuu.s us.i^h; KnuF>TH-iU';rH;r puMGF-*''ijT rA 
SUt.RUUfTiU-; Rt'«G7tALPHA2,AI.t'HA3,A) 

LU^ElVSitVJ Xf7,2(»0t), 1/(4,200!), 
ul '-'FAiSxrib l/CHGt-l ) vrihGl (4) ,FP(4) 

1U VK'.S104 crc4) ,CPi C4),CP2t4),CPi{4),rp3t4j 

uTmFGSIUU Vftp(7} rTi/(200l ) ,b'n(0) ,Z,!)C0) ,xn(«) « - 

'c A 2 / x7 A 37f ’ ’ ' "* “ 

Cf)riMf)N/Fi/ir:n,zn,HG 

C()irU',.iO/Fh 1 /PT, H 
COPMUH/Ay/EGU 
COf' .MufVRUNG21 /FP,VR 

UPTAL MF 
FEAt. 

A= A/5.0 

DO 200 Jl=l ,2000 
J1=2002-JI 
J=Jl 

on 55 1=1,4 
VPRf,! CI)=V{I.,J) 

cnr/nouF 

M=i 
i)=0-0 

DO 1 2 1=3,9 

zi(i)=za,u> 

COHTI54UE 
on 398 1=1,7 


Ito 


i 

2<»1 


44 4 


397 
2 3 3 3 


2 ? 7 . 3 : 

2 02 

2 v> 3 

20 o 

205 

7979 

5f> 

578 

399 
' 200 


aM(T) = X{ 1 ,vj3 

L>U 0=1,4 

V 0 H vj I ) •" V ( ' 1 n.i ) 

crjATifijif 


iM! 397 1=1,7 ^ i J 

/; r ^ ^ ^ ^ ^ ^ f 1 ' J - 1 5 ) 

III-, t lOuS 

{.01 23.33 1 = 3,9 

Z f I , J ) = H 0 . U ♦ t. X X- ( - Eij f 1 ) / ( F< G ♦ T J + z 0 c T ) ) 

<1=19 3 .9*( 1 .0 + U.a(n4*(T-4U,0l 
X2 = f.O 1 1,0 + 0. 00 1 4 ♦ ( T- 4 1 3 .0 ) ) / 3 uiM' . 
X3 = ( i 4.427 + 0,025 + CT-413,0) )/10nCi.0 
'*PP=(XC1 ,.T34X(2,J)+X{3,,1) + X(4,J3)/2.0 
vTij=XH'^FP + X2*K(5,JHX3*X(fe,J) 

!'f1 22222 T = 3,9 

TKC ( 1 .K0.4 J .UP. ( I,£0,9) IGfrrn 22222 

/,( l,‘n=Z(T,J)/VTU 

Cfli.TlMUf 



i)=U + 0,2 

CAI,!. F.uTh?(ALPmA2,40PHA3,J,T) 

Dtl 2''2 9=1 ,4 

C)'l ( -) = A + FF( M 

V(H0) = vn«GC]iJ-0.5*CPt (0) 

CHp r ir-iUF 

CAI,1, Iig'i9i2( ALPMA2,Ar,P4A3,J,TJ 

.in :?o3 t' = t , 4 

CP2f I'OsA + FPCtO 

VP(N3 = VnFGCfU-(t.b*CP2(0) 

CnOTlPlll' 

CALI. KgTN?lALPHA2,AT.PHA;3,7,T) 
iin 206 P=1 , 4 
C‘'’.3(0)=A+FP(!<i) 

VH{oJ=vrHG(i )-CP3CNl 
CU'.TINOF 

CAI-T. EUTf‘2(ALPHA2,Ar.PHA3, J,T) 

IM 295 M=t,4 
CiM(«)=A*FP(n) 

CP ( fi ) = (CP 1 ( P ) +2 . 0’t' 1 CP2 (N ) +CP3 ( 0 ) ) +CP4 f U 1 ) / 6 , 0 

V(^ tJlsVUPGCNl-CPrNl 

Cn.OTlNUF 



I p «»* 


.I'-uiii 1. 1.' Til ^f!i,vK An,jni?n’ kiouItfjJJs ‘ 
I’ti'jjT. r- F'^Til5(ALn,(A2, ALPHAS, J, F) 

i t r ^ / f ^ ^ ^ ^ ^ U , F P ( 4 } , V P 1 4 J 

> 1% ! i I • 1“ i/ A F C / j 

i'i{!%/Av/Lf:n 
'*irif‘J/TMiM<.i2t /l-'P, VR 


'iLATillu (,iF 


FXtLiT Pt-HIVATlVt;? OK 


■;AR(^) = l4.U**f2UPl-i729.0/T':4.042 + ALOu1.U 

2,';f ♦ ( H,0(>4ii)3-1757.85i/(T-S4.27 n j/ /A,..,o 
»'=i-T/C2,'.)*(VAPtF)-RT) ) 
v<Xl = P;4X2=P,*wX3=p;va4=R 
UAPt5}-PT) ) 

(Ui=v.>;GX2=y?GX 3=yjc;x4=..) 


ORi XI S-4.0 47C 3 ,n)*X( 1 t Jl + { 3,4 J /u.b) + ~X t 1 . 

l-Xf 2, J);X( 4, J) )*(;X t j )-ZX4 ,k?j-( j#(v( J)4Xf L a 

f JJ+2.0*Z(7.J)*X(2,JiT0Xl-if/,(7,J)/7.'; VI 
i»t )+‘vXl + Xf o. J) ) )-7(fi,J)*X(2,J) + { C7,(4,.T)/1.2'^ ) , p + ( ►•ro-A 

U . P-XC2..n-X(4,4))+WXt))-Z(3,J)*X(4,a) 

i)rix2=( (4.0*z(3,j)/0,5)^(-x{b,j) + (£r:D-xf Uj-s-xtZr-Ti-Af • vm. 
n)-2,0*Z($,J)*y(t f J)*GX2 + lC2.0»Z(7,a)}*CXi5,J) + A(2,J)*' a'/I)- 

1 C7(?,J)»X(1 ,J)/2.5F'^^2-Z(B,J)4X(i ,J) + ( F,.,o 

l + P-Gd-Xd ,.p-X(2,.J)-X(4,Ji ) 

an x3=4.l 

imx4=( f l.O=*=Z(3,J3/'*.5)*(-Xf5, J) + CEGI.-Xf l,J>-yc?,a)-Af i' 'I f * 
ic;x4 J )-?.(*i*7(5, J)*X(1 ,j)4Gx4+7.0^7,C7, J)yAf 2, J>*F:A4-.f4,f 7 ,.p ♦Yt 
l.I)/2,5 3*WX4f ( CZ(R,J)/1.25)»(-X{6,J) + fhr,ij-A( l,J)-XC?,4>-An,u 

HWX4))-Z(i,J3»X(l.J) 

DF2X1=Z(4 ,JUC (.2,0nC5,J))*(X(l,J1*r.Xl+XC5,J)3 )+2,0*Z(.7r.O*’<i 
.1 )-2,0*Zt7,J j*XC2,J)*GAt f({Z(7,J)/2.!>)»(Ar i ,J3»u.Al+XlH,a n j-,. 
1 J)*X(2,J)t( {Z(8.a)/1.25)*(-X(b,J) + (Lnu-Xf P J 

1WX1 )3-Ztq,J3/2,0 

DF2X2 = 2,0yZ(5,J)»X(! ,J)*LX2-((2,o*Z(7.dpy(XC2,vU»UX2 + A'b,>P > 
1 (Z(7,J)4X(1 ,J)/2,5)*WX2-Z(8,J)nCl ,J)P (7(B,vrj/j 
KFGU-xn , J1-Xt2, J)-X(4,J))*4X2))-Z(9,J)/2,C 
DF2X"3=0.0 , , , , 

l)F2X4=2.0*ZC5,J)4Xtl,J)*GX4-2,0*Z(?,a)»X(2,a)*GX4+tZ(7,a i»Af 1 
i/2.5)*WX4+( (Z(B,J)/i .25)*(-X{6,J) + fFna-Xf l,a>-Xl2,.P-A< 4,.p > 

aF 3 X 1 Sf r i ♦ Z ( 5 f J ) J * ( X ( 5 , J ) +X ( 1 , J ) ♦GX 1 J ) +2 . 0*7 ( 6 , .P ♦ A f 1 * • • ) 

UFiX2s2.U*Z(5,j5*XCl,.pXGX2 

DF3X3SO.O ■ ^ 

}>F3X4=2.0»Z(5,J)*Xtl,J)*GX4 

DF4Xl=-7(P,J)/2,0-ZC3,JPX(4,J) 

()F4X2=-Z(9.J3/2.0 
UF4X3 = 0,C» 

DF4X4s:-Z(9,J)/2.0-ZC3,J)*X(J,J) 

FP (1 ) = (2 . 04 ALPriA2*£G0* t ( EG0/X{ 1 , J ) ) -1 5 . 0 ) 3 / ( 22S . 0 ♦X ( I , J > * X » ! . 
l-(anxl4VUCn+DF2XdVR(2)+Of3Xl*7HC3HnF4xPVK(43) 

rP{2}=-‘2, 0'^ALPHA.3*)tC2f t 3 l 

FPn}i^(ijFlxUvn(t)+OF2X3*Vn(2}*OF3X3*V'R(3)^Of4XJ*vn^J)l 

FPf 43a-2 , f)*ALPHA3*XC4, J)"C0FIX4*VI3(1 )+i)F2X4*VP ( 2 > RUFiX^^'- R ‘ 3 > 

3+0F4X4>^VP(4)) 

HFtUHN 


a J + P-’.R 


:t 1 . T , 

*r A i 

1 V I ! . 

* 1 i im A ( 

' I -f itK / 
. I 

t> ,1. ) 


I I * 

I » Y t 1 . 

•* , I j 1 t 

I *'< 1 1 . 

I J - , 

, .p J* 

.on* 

1> , V) ) + 


Af 1 ,.P 
P »♦ 


j • 


1 1«J 


t- 

: 

¥ ■ 
r 


i'V' *1^'? 


IT 

;|y ' 

mV 

Ml 

mE: 

Ml 




. c 


)i)TlM 

K *3 j ( n i‘ 

.\n I'T' 
iSLn 
isw* 

■i3IT,i 

Glfl! 
jSiU 
NS in 


.C F~p'*snLVliiO”TriE’'"P"“nirFPH?rr Ty~“ 

.,-Knr'>'a NFTHOD ’ 

T F>inGl(ALPHA2, ALPHAS, A) 
Xf7,2O01),V(4,200t),VB('5j,AS(9J 
Zf9,20on ..iT2(101),HXUl-M J 
^ M4,4, K-i 1 ,P.SFG1 t4,4),PR(4, 4),pnHn(4,.U , 
i- CPU4 4),Cl>2t4,4).CP3{4,4) rpU4,4),(:'M t 


Tf 

i ) 


1UH/f.2/X/A3/Z/El/Ln,Z0,KG 
C'!M:1uM/An/TW/A?/y/AlB/HT2/A2r>/HTl 
CM CL)M/Ay/r'.GO/KF] /PT,H/B15/nKTP 
CMl 1iiH/An/P/Al4/PR/A21/FP/Al6/R 

A=a/(> T.iSATCuNL) 
iB'i 2'C:> Ji=l. ,100 . 

mi S3 [=1,4 
on S3 K=1 ,4 
nOPGlC L,K)=P( T,K,,T) 

S3 cn■^<T|UlJ0 

0=1 . 

0=0 

bB5S CCKf itiUK 

UfJ 55 1=1,4 
nn 55 K=JI,4 
FW(I,K)=P( X,h,J) 

P.irG(1,K)=F(1,K,J) 

55 Cnr'TlflUE 

L=C C J-1 )*2())-N + l 
Dfi 5959 1 = 1,7 
XtHT)=X(I ,L) 
iF(T,nT.4)GnTn 5959 
\[BCI)=V(I,L) 

5959 COin’INfiF 

T= lO'J ( L )-[>=>' C Tiy ( L) -T.1 (L-l ) ) 
iiClla)=56?0»EXP(-EU(l)/CKC='''T)+Z0(T) ) 

CDwTinur 
nri 44461 = 1,4 

k(I,L)=X(I,L3-a*CX(I,L)-XCI,L-l)) 
IF(I.GT,4)G0Tn 4446 , 

V( T,L)=V(I,L)"B»CVCI,L)-VCI,L-l>> 

4446 CnoTtNUF 

XI = 1 91 ,5* ( i ,0 + 0,0014* CT“4t 3,op /lOOO.O 

11 n ^ T •• 'll Q 

IFCCI.ioUl.OR.d.EQMDGOfO 2 

i CONTINUE 

B=a+ 1,0 /FLOAT CMNM,,) 


B 


1 ty 


IIU 

B ’ ' ' 

' 

203 
1 1 1 2 


1 U 3 


205 


i 777 
1114 


9735 


3U 


56 

200 

779 


C. 

c 


i)h 7<i2 K=ll4 

CM ( t , r ) rj \ r ,h ) 

,) = (. r , K J 5 * CP 1 ( i , K ) 

u I 20 3 i.sl *4 
COiC 1 , K ) 

^ .1 i 5 / = t’f n , F' ) -0 . 5 P2 ( I , K ) 

CM.c'; I'vr '.nAbr>!lR2,ALPHA3,J,T,,T3 

n J <00 • = I , 4 ■ 

(jp 2 00 Hst ,4 

CP3(X,**:) = A*FP(T.,k) 

tOUI,K} = pn,Xn(T,K)-CP3(l,K) 

Lll'iT 1 lot 

CALO FCT.M3tAbPjlA2,Al.PHA3,.l,L,TJ 

i>f.> 2')5 ls:1 , 4 

;Mi 2' 5 Kst.l 

CP4C J ,«s)=AiFp( I,K) 

cnoTIO’UF 
ufi 777 7 = 1,4 
UO 777 l..= t ,4 
Fl’F=:P(l ,1 ,*11 

lFCr>FP.nT. I .UF + OBMJOTn 779 

CtOiTIMiiF 
i 1 = f 1 + 1 

on 0795 1=1,7 
X{l, 0 )=XD(n 
lFi!,GT.4>GnTn 9995 
V(l,b)=Vli(l) 

COOtIMUF 

lF(M.b£.HtonGOTO 5555 

u = fj + l 
b=0,0 
0=1 

l.F(fl.LE,l9JGUT0 5555 ■ 

POHMATClX, M=' ,T3/1X, 'P(I,K,d)s',l6ri5.8> 
UO 5t> 1=1,4 
on 56 K=1 ,4 
FCJ ,K,.J-n=P(T,K,J) 

PC I ,K,d3=pnRGl CI,K) 

cnoTiNUF 

COfjTXNUF 

A=A4FL0AT(i7NM3 

RETURN 

f . Nl ) ■■ 

SUliROUfI}j“S‘'SDi;Vfe: ■'“P*’ Dl. ... , , 

SilHPUtlTl OF FQTNI ( ALPHA2 . AGPHA3 ^Jrl ,f) 

oTHFrision xn,2po|Lycl.20ou 
niMprisioM zN,2oot j reg(9|#zn(9) 

DIMENSION PRC4,4),FPC4j4) 
niUENSlON RCl,4£l01) 
l„) I M F r j SION ri T 2 C I P 1 ) , ^ , 

DIMENSION rX(4,4).XFC4 
Olf^ENSlDN SDHXC4,4) 



TP(1,4,101 > 
4,4) 


120 


C 'I / ^ ‘ / .A 2 / X / H 3 / 7 / C I / KOT^uTi^"' — 

( ' i )' ! if I i / A / / ^ 

CM", ni'i/ AV/E^fl 
C'' ^/.Al4/Pk 

CM'; i'ji</Ajf./P/nt7/PK 'fRR 
CM I ''V;/Alf>/»'‘J?/A19/FX,Xf=' 

CM j t/A7(VCL'HX/A21/FP/FFl/Pr,H 
C'm; nj’VA2f./rT1/Hl'i/nFTP 
CAM. PICCAr(J,L,l 1 

'-’Abl.. HAnLT2(J,l.,X) 

■HI 3 2A 11=1 ,4 

t'HP( 11 , 1 )=:TUR( 11 ,1 )*HT2{J) 

MHf.Tli.'-iP- 

ini,'?Mi'olCATTtJN MF I'^K mATPICFS "HR" AHf) 



III! 

3 

2ft 

K 

2 = 

1,4 


lift 

o 

'2S 

T, 

2= 

1,4 


H4i 

>l 

r. 2 

,E 

2) 

= TP!Hf 2 

)25 

* 

^ CO' 

M’ 

IN 

UF 




'""""c aI 

X 


XF 

Tj 



„ i|Hr 

A 

IP 

LTCA 

TTOM HF 

!188 

: !M) 

318 

v1 

2 = 

1,4 


lift K2=:l.4 

HFX(J2,K2)=0.0 
MM U6 1,2 = 1, 4 

MFX( J?,K2)aPFXCJ2,K2)+FXlJ2,L2)*PRCI.2,K?J 
3 to Cnt.TlI^UF 

318 Cni.TINUF 

C ' “ul'IplTc AT T ni~nr”"f }•!£”« ALICES •'XF"“u “j 

1189 .)fj 4 1R J2 = l,4 

DM 41 B K2=l ,4 
PXK(J2,K2)=0.fi 
on 41ft L2=l,4 

PXF(J2.K2)=PXF(J2,ft2)+PR(J2,L2)»XF(li2,K2) 
416 COwTInOE 

418 CncTlHUF 


CAM. HAFLTXTaLPHA2, ALPHAS, J,L,T) 

DM 4V6 1=1,4 
l)(! 496 K=1 ,4 

FP{T,K)=-Sl^>HX(l,K)-PXFCI,K)-PFX(I,K)+RMRiT,K3 

COrcriftUF 

HETUHM 

END 

MM ^F n ulP lI^^THl^SOLUrflft ^or'THi*^®'p " ^ aTR I x^ 

D t CFM S I ON X C 7 , 2 OU i ) , V ( 4 , 200 U , p C R , 1 0 1 ) , XH ( 9 ) 

DlfJiFMSinL •2(9,2001 ) ,E0(9) ,20(9) 

lilMiSsinii Srili(i)|DF2TX(4),0mX(4),fF4|X(4) 
DiHENSIDiv rTXa,4>,D0Cl,4),IOF'rp(l,|,lOl).HTiUOl) 
DTMENSlOAl FTP(4,43,RHCl,4)#PCl,4,i0l ),IRRC4,1) 
UlfAKHSION VAPC7) . ■ 


fY ^ K t A 



C'l //V/F.1 /FU,vn,RG 

C'j'r-inwAv/EClj 

Cn •' 'n; /A 14 /PH 

A 16 /R/A 17 /Ph,TPR 
f.:ij I' *■ /f* f '1 /t’X (,/ji 5/f)p-fp/^26/HTl 

V A 1 * ( 5 ) - 1 u , u+ ♦ f 2 1 , b 1 “ ^729 , C/T“ 4 , 042 ^ AliOGl <* f f T I / 7 fiu (I 

■iXi =p; wX^=[>;i*XJ=PrwX 4 =P 

r/( 2 , r* t v/\p ( s ) -pT) ) 

'»Xi= -,'s,X2=:y?GX3 = ..'rv,X4 = v 
;=\MPfS)=f f ( 37 29.0/1 i’+Tjj-4,042/T) 
a=vAP(o)M l7(>7.A53/(T-i 3.274)* *?j 
C = -(X(S,i,)*A + Xlf,,I,)*|j) 
u r=:(V( V APCS)-PT) 

J r=:CVC VApr 6 )-PT) 
p--r pT*A)/(2.0*((VAPt5)-PT)**?)3 
t:,=-(prtH)/( 2 .f*('iv/iPt 6 )-pT)** 2 ) ) 
six r'=:n;GA?’r=U?GX3T=n, *0X41-0 

>iXJ T=:l*;,? WX 2 T=H:? riX 3 TsEH«ix 4 T=E , v * 

■xnTRTTf+TTTp+oTOTT^TFTTfro TTTToi^To ^ 

;\2=Ou.n*t t .5 + 0,001 i*(T-4l3.5))/i00P.O 
X3r{ 19.4 22 + 0,O2*i*fT-4li.O)}/10U0.0 
Mps( X( 1 , 0 ) +Xt 2 ,L}+A ( 3 , 0 ) +X{ 4 , 1.3 )/ 2 .{; 

/!’=(( (iiPJ* I ,l>*O.0t)14+X(5,O)*t»o,f 1 + a(o,o) *0.020)'' » 

l + X 2 *(rf + X 3 *WT ■ 

i/r.i=:Xl*Mp + X2*X(5,0}+X3*XC6,L) 

00 t 1=3,9 

/PiCT)=bO.O*F:XP(-eO(T)/(RG*T)+ZU(I)) 

7.P ( I,. J) = r':oiT)*ZrMT 3 /(BG*T*T) 

lP((t.h0,4),()P,(I.P0.9))GnfD J 

;,PCT ,.I)s:i7P( I ,J)/VTn)-.( (ZN(I)/(V'fO*Vtf() )*Vp^ 

cnoTirniE 

oPiTXd J = (" 4 , 0 *x(t , 0 )-x( 4 ,ri))*ZP( 3 ,.n-c 4 .o/o,'=i)*( 7 t 3 , n*uT,.xi 9 . 

U.)*/.P(.3,a) ) + C 1,0/0.5)*CEGD-X(1 ,u)-X(2,ri)-\(4,L) )*(Z(3,3 •♦*.X iTt 

IGXI *ZP( J,J))-ZP(4, J)-2.{)*XCl,l.)*(Z(t),J)*GX!T+GXl*ZPt'>,vl O 
UZ( 5 .J)*CT+X( 5 ,L)*ZPd,a))- 4 , 04 Xn ,L)* 7 ,P( 6 ,J)+ 7 ,O*X( 2,0 7 i 7 , T 

i)*GXlT + GXi* 2 PC 7 ,J))-dU, 0 )/?. 5 )*(Z( 7 r 7 >*«d+' 3 Xt*ZP{ »-i 7 i 7 

l,J)*''iT+Xfp,L)*ZPC 7 /J))/ 2 , 5 -X( 2 fL)*ZP( 8 .J)-(Z(R,J>*V;r+Xl*..<;)* 7 »* 
103 , 0 ) )/ 1 . 2 ^+( CEGn-Xtl ,L)-XC 2 ,L)-X( 4 ,I..) )/1 . 25 .)*(Z{B,J) *» a) '» + 

aFfTX(L ’3 = r- 4 !o/O.S)*(Z( 3 ,wT)*GT+X( 5 ,Td*Z°C 3 ,J)) + C 4 , 0 *(PG( -Xi ! .f i- 
lX( 2 ,iJ«*xU,L))/ 0 , 5 )*CZC 3 ,J)*GX 2 T+GX 2 *ZP( 3 ,J))- 2 . 9 *XC! i 

1) ♦GX2T+GX2^ZP('3,J))+2.0*‘Xp,I,)*(ZC7,J)»GX2T+GX2*7P(/,J) * + 2 u+t7 
l( 7 .J)*GTix( 5 ,L)*ZH 7 ,Jd-bu,t,)/ 2 j)*(Sf 7 ,a)*’ 2 X?f + ..y 2 *>.PtJ.in 

l 2 ^f,)-X( 4 t))d. 25 )*c£{ 8 ,J)»WX 2 T+WX 2 *ZP( 8 ,j)) 

!iFrrX( 4 )s(- 4 . 0 / 0 . 5 )*CZC 3 ,vT)*Gf+X( 5 ,fd»ZPC 3 ,J)) + ( 4 ,U*CEGP-X« t .! *- 
tXd,t)-XC 4 ,L))/ 0 . 5 )*Czd,J)*GX 4 X+ 4 x 4 *ZPn,J))- 2 . 5 *X(l,Gi*tZi 5 . 1 
13 *GX 4 T+GX 4 npbr nb 2 . 0 #XC 2 .r.)»(Zn,JJ^GX 4 T+GX 4 * 7 P< ?,J) •-< xn .L 

1 3/2 S)*t7( 7, jbb4T+liX4*ZP(7, J>)-»(Z(8fU»XT+XCfe,I J *ZFt^ 

lvt))-XCl ,l 4 )'«'ZP( 3 ,J) . 




?T,aTT+ 2 -o*» 7 i^ 



f "! *. 1 + 

(2,L.)-Xt4,L))/} .'ysivCZfti,. 


r *' <:) 

* v ; x 2 T 
, s ) ^ - (. 7 
+■ Z P C H 

. P ' f-f ' J.X 

'TXC i) 
'TXCn 
i-ti,;X4T 
’ ( H J ) 
, 2 b )♦( 

s”xTTT 

. S , h )* 

i.TX(2) 

rr < c X ) 

ITXC 1) 

ifxTTT 
rr X ( 2 ) 
iTxcn 
I 'I X ( 4 ) 


= 2 . nfX ( ) , b ) 
+ f ; x 7 f 7 PC 7 ,.) 
C7,kIJ*WX2T + 
,.n )/i , 75 +t 

V 47 P ( 7 ,, 1 ) 

= 0,0 

=7.04xci ,r,j 
+ {7X44ZP( 7,J 
^ WT + X ( o , b )4 
ZCH,J)*Wx4T 

= 04 xTT 7 t 7 ' 

ZPt5,J))44, 
=7.0*X(1,CJ 
= 0.0 

=?,04X(1,0) 

rZzpT ^ 77 T 72 ‘ 

=-2P(9,J)/2 
= 0.0 

=-ZPC0,JJ/2 


t ( P d } l9l2Xt p2 npc sTj ) T=27o7xT I7i 

nZ(5,J)*GX4T4GX447,Pfb,J))-2,tn(2,i 
)) + ( X ( l , L }/ 2 . 5 >» CZ ( 7 , u 5 *’“ X 4 r + .. x .** 4 ,P 
J))/l.25t( (EG7-Kt 1 ,1 j-X(X,(,'i-v 
+WX4»ZPC8,J))-7P(9,J>/7'4 

*(7(5,J)*GX2T4GX7»7P(5,J) ) 
*(7.(5,J)»aX4T+GX447P(5,J) J 

75^47UT*fFrT7jT'”"^”’’^"”'””''”'~'" 

,0 

,0-XCl ,L)*Z.PC3,J> 



12 J 


1113 


: 14 


matrix li^Tn "TPP*' i'SATRIX 


or THF "HR" 

M.s:l ■ 

I'f) -iPS 1 1 = 1,4 
CHpri I ,K i -(sRhCKI ,T1) 

CRTf'itlF 
HKTIIHR 
!vll- 

‘TiiTs |~Hl<nTPfxNir’cALCULATes”THF”DeRrvAfT?es’’nF^F'‘ 4 

i;}iF''sir!S ?;f’;2"J{j'P*'2"<»),2PO,iou,2Bc<, 

uFxr 1.4 ) ,vap{7i 

oToF'SSJOM HTIC l6t ) ,OFTF(i , 4 , 1 0 1 ) , F0(9 ) , Z.Jf 9 ) 


H " f " 


cni,M.!iM/A7/v/t:i/fo,zn,R(; 

COMMll-v/AR/rnf) 

Cn;i 10>.|/A 1 ?/nFT/A26/!{Tl/FFl/PT,H/B15/L>FTP 

CAi,ri)i,HTTniM OF thf derivatives ok "f" w,h.r."T" 

HKAij '.HP . . 

VA1>(‘>) = 1 t).0*»pl ,6l-3729.0/T-4,042»ALO<;ior f) j//6o,o 
VApCf.) = in.u**(B.')641,ai-if57.8S3/(T-33.274>)/76().0 
A=;VAP(b)4C (37 29.0/(Tl'Tn-4,042/T) 

H=VAPfb)4( 1757.R53/(T-33.2?4)**2) 

(.’ = «(xib,L)4A + XCA,T,)»8) 

•;t=c/( vAp(5)-PT) 

W r=tV ( VAPCbl-PT) - w V . 



Kls:(lR,4 229O,02b*(r-413,0))/l0y0,O 
HP*(X(f ,b)+X(2,L)+X(i,b)+Xt4.LJ)/2,0 

VP=( ( (MP)* 191. 5*0,0014+X(5,L) *60, 640. 001 4+X(o,Ii) ♦0,025 • Oor.f^ » 

1+X2*GT+X3*WT 

V'T()=X1*MP+X7*X(5,TJ+X3*XC6,!:,) 

00 1 133,9 . 

ZfUI j5:bO.O*EXP(-EOtI3/CRG*T)tEO{n3 
ZP(T,vJ)3En{T)*ZNCI)/CRG*T*T) ^ 

ins 

CnMTiMUF ■ 

- -■* - - ■ZPt 6 ,J) + 2 .' 

,)-Xr4,u) 

,J)*Gt+X(5,0)*ZPlD,M) » 

1+2, 

1 I 

1 

iJPi j, I j ^ **»3 

U?l(K4)icCEr,O-XCt,li)*kt2,G)-XC4,L))/2*0)*ZPtR,JJ-X<4,i.l»Xit.Li 
1*ZPC3,J) ' •.-^..4,^ 4;;.^ . 

on 3 333 f ‘4 ^ 

DFl’P ( 1 . 1 . J ) «DPT< i f I t ' 



* s' " 
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ifss 


cn, T i ■'Hie: 

vi 1 ) ■ I SI t , U . , 

UH 1 -4 1 = 1 , J 

= 3 'i '‘4 v c f ,li) >!)(•?( 1 , 1 ) 

C'lrvr't!^'" 

,iTl f.JIsSllM 

-'i f A ".' v't , K , 1 , "T" 

;3i)npritrrL'.iF FftMLT 2 (J,L,T) 

[ST SinN !iT2f 101),Ff2(4),ZPU9) 

I'Sl.lFii.SiriN VAP(7) , . 

■cnnnu”7A27x7Al7z 

c:OM'i'iFv/A 7 /'//h:i/po,zn,RG 

CfJMMClN/Ay/KGU 

Ci)f.M(iy/A IH/H't2/FF1/PT,H 

CALCtHiATinii UF TIIF SfiCuMU UFHI V A I* I VRS vjF "F'* W.R.^’,"T” 

Kl'M, FP 

^AP(‘j) = l 7 . 0 *T ( 21 , 61 - 3729 , 0 /T- 4 . 042 FfiLnG 10 fT) 1 / 7 m ).0 
i/AP (6 ) = t 0 ,O*t(H,O 64 l 03 -l 757 ,R 5 i/{T- 33,2743 J/ 76 U.I* 
A=VAP(bU?( 3729 , 0 /(T)^tJ)- 1 - 042 /Tl 
H=VAP( 9 )-* f l 757 ,R'i 3 /(T- 33 . 274 )*» 2 ) 

C = -(X( 5 ,L)*AtXC 6 ,T.)* 6 ) 

HTsC/C VAP( 5 )-PT) 

.•vT=C/(VAl>(&)-PT) 

G=VAPC 5 ) 4 ( (- 2 , 0 * 3729 . 0 / CT*<' 3 }) + ( 4 . 042 /(T» 42 ) O + f l 3779 .o/ 
1 - 14 , 042 /T))*A 

|J='/APC 6 ) *(- 2 . 0*1 757 . 853 /C(T- 33.2741 4*3 i ) + ( l 7 S 7 .P 57 /f {T-i 

i**?) )*B 

GT 2 =(-X( 5 ,L)*G-X( 6 fl.)*W-^T*B- 2 . 0 *GT*A)/(VAPfS>-PTT 
'VT 2 =(-X(!>,G)*G-X( 6 ,L)*!«-GT*A- 2 . 0 **lT*p)/(s/AP(o)-Pf ) 

TnT9T:577T:^f3r(^TPTF?TOTi7rs7^^ ^ ^ 

1 +X 2 *GT+X 3 * 1 «JT , 

VT 0 =X 1 *MP+X 2 *X( 5 ,L)+X 3 *X( 6 ,L) 

00 t' I **** 3'' 

zfu r 3 =60 1 o*EXP ( -EG a ) / ( Ra*|Hza( u ) 

ZP(I,J 3 =l?: 0 (T)*ZN(I)/(gG*T*T) 

IF((t,E 0 , 43 , 0 P.CI.Eg, 93 )Gg |0 1 

ZPC I » J) = ? 7 ,P?I r J 3 /VTd)-((ZM( J)/(VT 0 *yT 0 ))*VP) 

cfntinuf 

GlIT'T 2 ' '■ '"■ ' ' 


t ' X, f , I ' 


4 I 


i . / 7 '4 I 


I 0irO *0 i 


inn 
1 C) 



<JXls;.j;GX2sp07GX35'yrGX4«a 

FKT7TTs-TZ7i^|n7?IIlIlT^ 

1 udrCyJwxitxis.tifljZJs 

pill 

l + iFGO-XC 1 ,L)rXC2-fl‘l-X;(4,, 

FX(3#1)»0,0 

rx ( 4 ; t ) * t ? 4 . Ot'lSC 3 , J:5-/»*f ) 


1 ' Mf ) I’l-y j, j j j-Ai ’ 

7,:.,p: sW z' ^ 

1 (" '> f 7 * i ^ / * v ? r ^ *v.V^ ^^■*'2,0^X(l|L3^XflfG)^Hp(<5»iJ)“’2,P + Xi 

t (. / ,( / j. iH ♦M 7 + X ( "S , G) ♦HP (7 ,tl }+2,Q*Q’f *7,9(7 iJ)) + (X(1 
iJ)*..’!> + Uo,bUHF{7,J)+7 0rr|?ZP(7,i)>-xt7!l. »A 

Fr2?4l = r fFGO-Xtt ,i)-XC2fL)-XC4,l.))»HP<9,J) j/2.0-Ar4,b'> ♦> ' 

1 ( 3 , vO . 


.S')n = 

DO 3 
sn.ir 
Cn;’T 
'iT2( 
tlT'^f 
({(•’Tfl 

IvJi) 

mmti'tmfiiimmm 

miG 

Siii)K 

UlIlF! 

ni'-tp; 

iMMK 

uX'^f; 


3 / 3 Ts! t 4 
= si!M4 v(i ,LnFT2CI J 

V Li'JijF 

( J)=F.liM - 

f.Osl ,0/HT2Cil) 

IH'-J 

rsuHRnTifTNrcA"CT!:rfEs*i^rf>mvAfTiFns’of”^^ 

•bJUTlbE FXF(J.L,T) 

^HSION K(7,206i KvC4r2001 j 
'bM.Sinu ZC9,2001 ) 

*:HSini. FX(4.4J,XF(4,4) 

VAP(b 


COMffJH/A^/X/AS/Z 
Ci!.i.''lOM/A7 2V 
CnilMON/Ay/EGf) 

CnMM()N/A19/FX,XF 

C(l»1MUN/FFl/PT,H ,,:. • ,, 

? A 1’ lol =I 0 * I ♦ ♦ ( IT sf 1 1 o|^?f ^ ,851/ f 1-11 . 1^4 ) } 

F=FT/(2.a* (VAP(fe)-Pf)) 

'.'/X 1 = P ; WX 2sF ? WX 3sp y WX4=P 
GsFT/ta.O^CVAPCSX-PT)) 


nGT0(T>i/76i>.t;> 

.274)}/7f;.),0 


TOtfTTPoTfTTTivT^TrTTrG-'-^ , j ^ 

M(Kon(5:.J)3*txh;T.J + XU ,i.74 

7lj)fXC2»ti5*GXl-C|ZnrJ3^^.‘’ )♦ 
2JlH((2C8,vJ)/ 1.25>4(-XC‘ ' + i 

iM^«Ga-XCl,G7-XC2,I.-)-Xt4,h-> !♦ 


ooo o 000 
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c. 

c 


301 


10 

21 

7.7 

n 

24 

It 


,1,! 
1 I } 

r I 


1 + . K4)5-Zf J^jU^U “ 

I _* 5 X J J *** L ' f t,l / / in V 

Fx f ;• , 2 ) = 2, ')t /, ( 5 j.j ) *K ( 1 , L) *Gy2-( (2,o*2 ,( 7 fill ) ♦ (. V I •;> ,r, ) » ,x^ t * < s , (. 

1 I f ^ ' V f 1 i ^ ' f ^ ^ vf 4 i ^ ^ ^ + t f /i f "I , J ) / I . 2‘J ) i ~X > n . 

1 + { ;.«Ti i-X ( 1 , 1 ) - ^ f 2 , U t -X 14 ,T. J i ♦ KX? ) )-7. (9 , .1 j /2 . " 

FXC 3, U = 

rX(4, ? J = 2.<>=^/:fS,J)*X(1 .L)t3X4 

FXt t , l)=-Z('^,.7)/2.<!-Z(3,J)tX(d ,t.) 

FX(2,4J=-Z(9,.7J/2,0 

FXC.}, It = ’).!) 

f X( 4, i) = -Z(9, J)/2.0-’Z(i,J)»y(l ,r.j 
FXf, 3,2) = 0,') 

FXti,2js:2, i*Z(bjJ)*y(l ,LO*Gy4-2,04'ZC7, J )*X( ?,T. >*Gy4+(.7t / ..f i*a( 

I I ,6)/2.f>7*l.x7 4-( C7CB,J)/l.25l*t-X(h,L,)f Cc<lU-A(i ,Lt-XC2,F )-Af^, 

.1 1 J*iJX4 ) )-Z(‘J,U)/2.0 . . 


■•^ATRtX 


rK.’\:i;3fniSt (jf thr 
JM 301 Kir 1,4 
i>n ^ui 11 = 1,4 
XRf. ri rKl l=FXtKl , U ) 

COM'tfaiR 

KRrnPK) 

RfiP , *• 

TFf I S^SMBFCIIT 1 6R“L?ULlT?S™K'"s?CUHu”DRKT7nT?Rr"uF ^ V ^ 1 , • 

TflNTA^ W.R.T,"x'' 

SnRRUU'riwF HAMLTXC A6t>HA2.ALPHA3,3,U,‘n 
OI".i:i-SinM X(?, 20ul 5 ,VC4/20»Jt ) 

UT’.-RGSin,. Z(9jr200l),Sr>jXC4,4J 

i)Ti'F.^S|nK flX2C4,4),F2X2C4,43,F3X2C4,4),PMK?t.4,4l 
DlFFXSiex FX2(4,4),Sl)HX(4,4) 

OT MR?'*. SI On V4PC73 «■ V 


nr ili.T i.t‘ j fl 


.0 


’c n" m^n 7 A 27 X / If Tz / a T7?7I?7RGn' 

CnHMUM/A20/SUHX/FFl/PT,H 

’CAlc!lLAfIniraF"f?P:"'flco6crfJfHf7ITT7f?r"o?'”'’^j^TI'’7^I’ 

W,R.T.'*X” 

00 10 1=1 ,4 
00 10 K=l,4 
SnjxCI,K)aO. 
cni.Ti’^iiF 
UD t 1 1=1 ,4 ^ 

siMX(ljT5 = (^!odLPHA2tEOOF{3.0^RGO/X(lVLj-2.0*l5.0j>/(22h.u*A(» 
U,)*X(I ,L1#xh,T.3)FSDJXfIrI) 

GOTO 11 

firiox(I,I) = 2,0*ALPHA3 

goto 11 

SDjX<t,T)=0,0 

goto 11 

3fMX ( I , I ) *2 . 0*AhPHR3 
COMTIftOE ». ^ 
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I SitTUf-.’U HFRT 

':i )s:1. *,0*4 c^l .6 1-37 24.0 /T-4, 0424 u ( r ) j / / <S'j . ' 

^ Ap C b ) = 1 . I,! ♦4 ( H , 06 4103-1 75 7 . 85 3 / (T«33 , 27-1 ) j / 7ft(' 0 •' 

P=PT/{2,0llVMi’(03-PT)) -ijjfir,,, 

OX t=P? .7/2=;PrwX3=Fy-a4=p 
■) = PT/( 2,0*1, V/^P(5J-PT) ) 

<;x J =v7 '’Xl>=:0;f;x.3=:0;GX'l=Q 

’™t^””77j7=~:f7|i?^jf-jji^jj72Y67jT-?47o*z7I7oT7o7s7T7’77( 

I '-xt - C { 2,047(7 ,J)/’2.5) + (2.U*2(8, J)/1 .25) )*».'X1 

n + C 2 , '.j + 4 ( h . J ) ) 1 *f:x2- C t Z 1 7 , J) /2 . 5 ) f c 4( 8 , d V I , 25 ) 7 X 
I C « , d ) / I , 2 5 ) ♦ V-VJ 
F1X2( 1 , 3 7 = 0,0 

Ji;;v4-{X(rj,d)/1 ,25)*WXt-C(Z(7,d)22.5U(4(8lj1/l ,25) 7^ A. 
FlX2U,nsFlX?( J ,27 

FJX2(2,? j3(.(-4,04Z(3,d)/0,5)+4,0*Z(7,d)l*nx2-(2.')*Z(.S,. 
1.U.X2 

F1X2(2, 3)=0.a 

FU2C2,4}si-2,»j4ZC3,J)/0.5j*GX2-C(2.0*Z(3,d)/0.5 7-2.U*: 
l0X4-(.Z(8,d)/!.25)4WX2-Ct{H,d)/l .2554WX4 
-PI 20 1=1,4 
K1X2C3,T) = ( 7.0 

continue 

FlX2(4,t ;) = F1 X2(1 ,4) 

FlX2C4,2j=f- 1X2(2,47 
K2(4, 37 = 0,0 

FlX2(4,47=(-4,0*Z(J,J)/0,5j*GX4-(2.O*ZC8,Jj/1.257*wy I 

'? 2x 2 TT7TT= 47 o7z( 57 jTtGxT+4 7o7z c 6 7 jT+727o ♦ flTTd T'27 5 ) *' ‘ 


I, 71^ Ji/i 

F?l2ri!2) = 2,5 *Z( 5,J>*GX2-2.0 tzaf 7)*GXl + UZf7 ,dl/ 2.5 7- 

II, 25))twx2-ZfS,d)-(ZC8,J)/t,?5)*^Xt 

F2X2n ',4j = -ilo»Z(5,a)»GX4 + (CZC7,.I)/2.5)-(Zf d,d)/l.25) 7* 

U/1,25)1'WX1 

F2X2(2,l)*F2X2(i,2) 

F2X2aj2)=-4,oiz{7,J)*GX2-(2,0«=Z(8,J)/i.25)*ftX2 

F2X2(2,37=0,5 

f2X2(2,4j=-§.04=Zt7,J)*GX4-(Z{8,J)/l,25j*NX2-(7(8,J)Ai,: 

170 30 1=1,4 

F2X2(3,I)»0,0 

CONTINOF 

F2X2(4,1 )=F2X2Ct ,4) 

F2X2{4,2)=F2X2(2,4) 

F2X2(4l47cC§!o»Z(9,J)/1.25)’^iilX4 _ 

■mTrT7TT5 4:PzT?73Tff|T7¥:TO(T731 

F3X2(l,2Ja2.0*Z(5,J)*GX2 

r3X2(l,3)=0.0 

F3X2hl4) = 2.0*ZC5.J)*GX4 ' 

F3X2(2,I )»F3X2C1,2)„ n 

FiX2a,2 3 = a,0jr3X2(2,3)»0,0?r3X2(2,4)=O,0 

00 32 1-1,4 . 

F3X2C3,I)*0,0 

COHTlNOE 

F3x2C4^t )»r3X2(l,4), 


1/2.5 7 ■ 


..T< / 


2 1’ ’I ’ 



12 b 


14 


■14 

15 

4ft 

42 

49 


50 

C 


iBA'n 1 ,2 J=»-3X'>C2, 1) 
^B'<2C 5 » B^sO.O 



.V^ 'I 3 4 f = l,4 
■)<' 3 4 h=U4 
CKi r I ,ko,i ) 

ri.!io,4) .AOD, 
!-’4X2(. i:,K jS'J.O 
CiO.l'I'JijK 

,4) = -7(1,J) 
i'’4x2(4,l J3:F4X?(1 ,4) 


tK,F:o,4))GOTO 34 

))CjTa 34 


the adjoint VAHTA8UF3''Tu"‘ff?r 

A 




-IJi.TIPblCA riUN OF 
VTb OF "F” F'.P.T, 

Of 42 1=1,4 
O'-i 42 1*1,4 
JO 42 K=l,4 
';nTnC4l,44,4ij,46),4 
FU2(T.!<J = V(M,jL.)%FlK2Cl,K) 
goto 42 

f 2 X 2 (Ie?<l = '/(M»L)'t'F 2 X 2 (I,K) 

GOTO 42 

F1X2CI,M = V(M,1.)*F3X2C1,K) 

OOTO 42 

F4X2( r ,K)=:V(M,b)*F4X2(I,R) 

goto 12 

COoTlNlJfi’ 

DO 19 1=1,4 
on 49 f=1,4 

FX2( I ,K)=F1X2(1,K)+F2X2(I,K)+F3X2CI,K)+F4X2(I,K) 

COMTl^niF 
DO 50 1=1,4 
DO 50 K=l,4 

Sl)HX(t,K)=snjXlI,K)+FX 2 Cl,K) 

coi'n’iNuF 

I'YPFJ^J, C(SPHX(T,K),I=1,4),K=1 ,4) 

FO'l) ' % . ■ 

T HI s*^ u bTo uf srJ b?l s"*T Hi r^FF R F jtII u’TouIr 1 ni s“ uF r F u'* K IT 

-E-KUTTA METHOD 

iiflHHUUTiNE RyNG4(ADPHA2, ALPHAS, 5) 

DIMEDSION XC7,200l),V(4,20QU 
DlMENSIOii ZC9,200U,TN(200U 

DIMEHSBIN U(4,1,101),OORG1 C4,U,vR14,13ignHG(4, 1) 

UTflEHSIOH C0U4. U ,c62(4. ' ' - - -- 

DIOENSION C0C4,1 ),Fgf4,l ) 

l)lM|flsfoL H?2c|6f IBITiClOl ),0FtP(t,4,tOt ) ,K(1,4, 101 ) 


(4, 1),C03{4,1),C04(4,1 1 


’c 0 M H DnTA 27x7 A 3 / Z” 
C0MHUN/A7/V/A6/Tf 



COMMON/A26/HTJ/A18/ 


6/R 


D0^30S‘^ai»tf 100 

' ‘ ' -.v 

■ -uv ..r. 
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.!=ul 
t = ! 

.;n s'i h=i ,4 

.''“.'CU (r , I )=CIK,1 ,.T3 

LnM'i'aU*'. 


GHn 

!.M) K=1 , 4 

V^'(r r n=U(F,I,J) 
jnr<r.(r.,i)=y(K,i,j) 

Cn.'.iTl^iUF 
fj=:( (J~l 

rsTfj(L)-u*(TN(L)-T?J(L-I) ) 

TaTWU.) 

i,»o I r=:3,q 

cnl r c -Fo ( n / ( rg^t) + zo( r ) ) 

xisi')l ^'i'K I .0+0,0014*(T»4t3,0) )/1000,0 

X2 = f;'i.bf{ I .f5 + O,5014*(T-4n,5))/lo00.0 

U.*;* + 3.0 )J/1000,0 

"=(X(i ,L) +X(2,l;)+X(3.ri)+X(4,L))/2.0 
\/Tri=:Xl’*'MP + ;(?*X(5,L)+X34XC6,L5 
un 111 1 1=3,9 

;IFC(T h0.4J.iJP.(I.f;0,9)>GnTO Ull 
[i=B4U,2 

CAl.L EyTW4(Ai4PHA2,ALPHA3,J,L,T) 

l»n 202 K = l ,4 

CUl CK, n = A*FufK, I) 

5HCK.i)=onRGCKj )-0,5*C0HKa) 

niiir.iNfj!' 

CALfi E(2TN4(ALPH A2, ALPHAS, 
un 203 K=1 ,4 
C02(K, C)=A»Fy(K,I) 

^)RCK,!)=QnHG^K,t J-O.S^COaCK, J.) 

CALL EQf<atALPHA2,ALPHA3,J,L,T) 
on 206 K=1 ,4 
C03(K,I)=A4‘FOCK, I) 

:JP{K,l) = 0nPG(K,l)-C03CK,I) 
cnwriouE 


CONTIOUE 

CALL Egii>i4CALPHA2, ALPHAS, J,L,T3 
on 205 K*l,4 
Cg4CK,l)=A»Fg(R,I) 

CO(K, I)!?(C91(Ka)+2.0*(Cg2CK,I)4-C03tK, t))+C04CK, n)/b.O 
9(K,X,J)=OORGCK,i3-COCK,I) 

continue 

lF(H,LE,5)(jnT0 5555 
■^■*0 + 1 
8 = 0,0 

IF(N. 1^6,19)6070 5555 
DO 5p K=l,4 


QCK,i,a-n«gc 5 

CONtlMUfi. ■ ' ' ' 


mkii::::- 




,• ' 4 i.’n * > 
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' C WWM** 

202 

ti 

9 

n 

10 

n 


12 

c 


15 

u 


16 


A = A <‘'>,0 

ttH rui-^ i 


i> u H r ^ n r [ J ? i K T u s 0 ^ii YFf E We W f T A ^ 


in.'KM.Slfi XC7,/!0Cl),y(4, 20011 

iirv;|.-„!.S[n|i Zf9,2001 J 

;jrjK-ir,in.i ooi4,i ),f.Q(4,1) 

'! I ’ ’T" r J pi! nH n ' f S ' S’U ^ J V n 1 f i , , 10 1 ) 

oT'il M.Sir-r' H'rs(4, 1 ) ,DFTn( 1 ,4) ,PHMfc;WC4, 1 ) 

i'l ‘ ^' hSlOf' XF(vi( 4, 1 ) , R(1 ,4, lul 5 ,RP( 1 # 4 5 ,1 RhC4, 1 1 

ntM-w.Slf'W f-X(-i,4),XFC4;4LEnC9),ZnC9) 

'cnr!inw 74 27 r 7 Al 7 z 7 Fl 7 io 77 :o 7 KG' 

(:n;r).)H/A7/v 

C nn MO W / A O /EGtj / A 1 5 /DF t / A2 6 / HT l / B 1 5 / DFTP 
CO-v^'1l)N/A16/R 
Cflf 1<,i"i/A18/HT2/A10/FX,XF 
Cnt.,’inM/A22/0fi/A2 3/Fa/FFl/PT,H 


■c^rFxPuTITfT 

Uis 202 T6al,4 

DFTtl , 15 )suFTP( 1 ,I 5 .J) 

COriTlNUF 
K = J 

1)0 ,H 1 = 1 , 4 
HFUK,T3=F(K,|,J) 

TRKCIrKlsRHCK,!) 

CO'vTINUfc: 

DO 0 1=1,4 . , 

HT5(I,K)=TRR(I,K)#HT2(J)*HTICJ) 

COhtlNUE 

nULTlPtTCATIQN OF THE MATRICES "DfT" AOD "OP** 

00 10 1=1,4 

DFTM(K,ll = HT2(J3=i‘DFT(K,l) 

COMTINUK 

nFTO=DFTj+DFTN(K,I)*QR(I,K) 

continue 

UO 12 1=1,4 

RRNEW C I , K ) sTHR C 1 1 K } *DFT0 

continue 

MUI.TIPLIjCATICjN or the matrices "FXF" and "y»« 

11=1 

on 14 Kl=1,4 
XFOCKl ,li)»0.0 

1)0' ' f !§' ft 1 S'"l *4' ' 

XFy(Kl,ln=XFaCKl,Ii)+XF(Kl,Ul)»aR(Ll,Ili 

CONTINUE 

CONTINUE 

FO(i£K)=ifl€I,K)+RRNE«CIfK)-XFQ(I,K) 

CONTINUE 

RETURN ..... , 

. .. _ 

lAL EOUA'‘I•)^ 


•ffll 

FORNA 

5 UBR 0 ' 



A 2 >Al»f»BA 3 fA 3 


», .f !.>..» ..h' .*- m...’’ . 

O’. *dv 


131 


■ jUK' 

Fi'is'in:. 

FuSlflM 

iF'JSTni-, 


Xr/,2001),V(4,5>0()1) 

Z(y, 2 (Hii) .Twt?ooi) 

r 1, n * f 4 f 1 4 \ r\ r ♦ 1 - * V # .4 


J.TWr/OOU 

['fJ*'5H‘ht5#DFLL.X(4,l,lOU 

Pr n ' } 1 1 f 1 < CD3 ( 4 , n , CD4 ( 4 J ) . 

nTuion,HrniGi j,iSF'fPu;4,io{j 


b:n(9 j ,zn(^) 

){Cnt 4 ,i X 

f ^ ) 


Cnf!M(i:;/A2/K/A3/?VRl/Fn7zu,R?5 

C:!);iMU;)/A?/tf/A6/TW 

cnh;*nii/A'4/fc;cii) 

C'V,v.rj;i/A2'-V!’Er.LX/A27/l)ET<l.Xi< 

CUtintli'M/A2R/TttV'-X/Fri/pT,H 

CnMMrjr)/Ai6/R/A24/0/A26/H'ri/Al8/Hl’2/pl5/DFTP 

^ ^ 

As:A/‘>,0 

I'H 100 J 1 = 1, 100 
vJ=U1 

i.i=l 

on 00 K5'1,4 

unfjGl (K,L)=nELLX(K,L,J) 

COMTtMUE 


C 0=0.0 

f’OHTLOUK 
DO 65Mv=l,4 

DELljXR(K,L)*l)EbI.X(K,U,J) 

DELXHG(K,r.) = S)EI.LX(K,L,J) 
tj5 CnoTINilE 

T=TWCM) 

z( j;iMaio!o^l:xp(-EQ(n/(KG<'T)+zo(r)) 

mil xt = 10lli!i*( l ,0+0,0014*(T-413.0))/1000,0 

X2=60,o*( 1.6+0. <loi4*(T::4i 3.3)} /1 000,0 
X3=a3, 422 + 0, 0254Cf-4i3. 0)1/1000.0 
i'1P=CX(I,M)+X?2,M)+X(3,'4)fp4,H))/2,0 
VTD=X1*MP+X2*XC5,M)+X3*XC6,M) 

XFC(iPo.4k6r.(I.EO, 9))GOTO 22222 
z(i,j)sZ(i,J)/vTa 
22222 CONTINUE 

C 1)=B + 0,2 

■ CAbU EyTN5CJ#M,T) 

DO 902 K=l,4 

CDlCK,L) = A*TqWX(K,L) , _ c*rnirK 

„ DEbriXR(K,L)=0ELXRo(K,b)+0,5’«'Cf)i(K,b) 

902 continue 

U§0bliflW*l)=8KlR6W,l)t<).5*CD2CK,I-) 

903 CONTINUE v. • .'«tO ■ ■ '• 

' CAUL r.aTN5(afM,T).':..A:;: 

DO 906 ,K»1,4 . < 














